diff --git a/python/README b/python/README new file mode 100644 index 0000000000..742d5e6f62 --- /dev/null +++ b/python/README @@ -0,0 +1,47 @@ +This directory contains Python code which wraps LAMMPS as a library +and allows the library interface to be invoked from a Python, either +from a script or interactively. + +Details on how to build and use this Python interface are given in +doc/Section_python.html. + +Basically you have to extend the Python on your box to include the +LAMMPS wrappers: + +python setup_serial.py build # for serial LAMMPS and Python +sudo python setup_serial.py install + +python setup.py build # for parallel LAMMPS and Python +sudo python setuppy install + +but there are several issues to be aware of, as discussed in the doc +pages. + +------------------------------------------------------------------- + +Once you have successfully built and tested the wrappers, you can run +the Python scripts in the examples sub-directory: + +trivial.py read/run a LAMMPS input script thru Python +demo.py invoke various LAMMPS library interface routines +simple.py mimic operation of couple/simple/simple.cpp in Python +gui.py GUI go/stop/temperature-slider to control LAMMPS +plot.py real-time temeperature plot with GnuPlot via Pizza.py +viz.py real-time viz from GL tool in Pizza.py +vizplotgui.py combination of viz.py and plot.py and gui.py + +Run them with the following input scripts and arguments: + +trivial.py in.trivial +demo.py +simple.py in.simple +gui.py in.gui 100 +plot.py in.plot 10 1000 thermo_temp +viz.py in.viz 100 5000 +vizplotgui.py in.viz 100 thermo_temp + +You can un-comment the Pypar calls if you want to run these in +parallel. + +Each script has more documentation at the top of the file that +explains how to use it. diff --git a/python/examples/demo.py b/python/examples/demo.py new file mode 100755 index 0000000000..be38d30167 --- /dev/null +++ b/python/examples/demo.py @@ -0,0 +1,66 @@ +#!/usr/local/bin/python -i +# preceeding line should have path for Python on your machine + +# demo.py +# Purpose: illustrate use of many library interface commands +# Syntax: demo.py +# uses in.demo as LAMMPS input script + +import sys + +# parse command line + +argv = sys.argv +if len(argv) != 1: + print "Syntax: demo.py" + sys.exit() + +me = 0 +# uncomment if running in parallel via Pypar +#import pypar +#me = pypar.rank() +#nprocs = pypar.size() + +from lammps import lammps +from lammps import LMPINT as INT +from lammps import LMPDOUBLE as DOUBLE +from lammps import LMPIPTR as IPTR +from lammps import LMPDPTR as DPTR +from lammps import LMPDPTRPTR as DPTRPTR + +lmp = lammps() + +# test out various library functions after running in.demo + +lmp.file("in.demo") + +if me == 0: print "\nPython output:" + +natoms = lmp.extract_global("natoms",DOUBLE) +mass = lmp.extract_atom("mass",DPTR) +x = lmp.extract_atom("x",DPTRPTR) +print "Natoms, mass, x[0][0] coord =",natoms,mass[1],x[0][0] + +temp = lmp.extract_compute("thermo_temp",0,0) +print "Temperature from compute =",temp + +eng = lmp.extract_variable("eng",None,0) +print "Energy from equal-style variable =",eng + +vy = lmp.extract_variable("vy","all",1) +print "Velocity component from atom-style variable =",vy[1] + +natoms = lmp.get_natoms() +print "Natoms from get_natoms =",natoms + +xc = lmp.get_coords() +print "Global coords from get_coords =",xc[0],xc[1],xc[31] + +xc[0] = xc[0] + 1.0 +lmp.put_coords(xc) + +print "Changed x[0][0] via put_coords =",x[0][0] + +# uncomment if running in parallel via Pypar +#print "Proc %d out of %d procs has" % (me,nprocs), lmp +#pypar.finalize() diff --git a/python/examples/gui.py b/python/examples/gui.py new file mode 100755 index 0000000000..4fa3c36d47 --- /dev/null +++ b/python/examples/gui.py @@ -0,0 +1,109 @@ +#!/usr/local/bin/python -i +# preceeding line should have path for Python on your machine + +# gui.py +# Purpose: control a continuously running LAMMPS simulation via a Tkinter GUI +# Syntax: gui.py in.lammps Nfreq +# in.lammps = LAMMPS input script +# Nfreq = query GUI every this many steps + +import sys,time + +# methods called by GUI + +def go(): + global runflag + runflag = 1 +def stop(): + global runflag + runflag = 0 +def settemp(value): + global temptarget + temptarget = slider.get() +def quit(): + global breakflag + breakflag = 1 + +# parse command line + +argv = sys.argv +if len(argv) != 3: + print "Syntax: gui.py in.lammps Nfreq" + sys.exit() + +infile = sys.argv[1] +nfreq = int(sys.argv[2]) + +me = 0 +# uncomment if running in parallel via Pypar +#import pypar +#me = pypar.rank() +#nprocs = pypar.size() + +from lammps import lammps +lmp = lammps() + +# run infile all at once +# assumed to have no run command in it + +lmp.file(infile) +lmp.command("thermo %d" % nfreq) + +# display GUI with go/stop/quit buttons and slider for temperature +# just proc 0 handles GUI + +breakflag = 0 +runflag = 0 +temptarget = 1.0 + +if me == 0: + from Tkinter import * + tkroot = Tk() + tkroot.withdraw() + root = Toplevel(tkroot) + root.title("LAMMPS GUI") + + frame = Frame(root) + Button(frame,text="Go",command=go).pack(side=LEFT) + Button(frame,text="Stop",command=stop).pack(side=LEFT) + slider = Scale(frame,from_=0.0,to=5.0,resolution=0.1, + orient=HORIZONTAL,label="Temperature") + slider.bind('',settemp) + slider.set(temptarget) + slider.pack(side=LEFT) + Button(frame,text="Quit",command=quit).pack(side=RIGHT) + frame.pack() + tkroot.update() + +# endless loop, checking status of GUI settings every Nfreq steps +# run with pre yes/no and post yes/no depending on go/stop status +# re-invoke fix langevin with new seed when temperature slider changes +# after re-invoke of fix langevin, run with pre yes + +running = 0 +temp = temptarget +seed = 12345 + +lmp.command("fix 2 all langevin %g %g 0.1 %d" % (temp,temp,seed)) + +while 1: + if me == 0: tkroot.update() + if temp != temptarget: + temp = temptarget + seed += me+1 + lmp.command("fix 2 all langevin %g %g 0.1 %d" % (temp,temp,seed)) + running = 0 + if runflag and running: + lmp.command("run %d pre no post no" % nfreq) + elif runflag and not running: + lmp.command("run %d pre yes post no" % nfreq) + elif not runflag and running: + lmp.command("run %d pre no post yes" % nfreq) + if breakflag: break + if runflag: running = 1 + else: running = 0 + time.sleep(0.01) + +# uncomment if running in parallel via Pypar +#print "Proc %d out of %d procs has" % (me,nprocs), lmp +#pypar.finalize() diff --git a/python/examples/in.demo b/python/examples/in.demo new file mode 100644 index 0000000000..8eb9020139 --- /dev/null +++ b/python/examples/in.demo @@ -0,0 +1,26 @@ +# 3d Lennard-Jones melt + +units lj +atom_style atomic +atom_modify map hash + +lattice fcc 0.8442 +region box block 0 10 0 10 0 10 +create_box 1 box +create_atoms 1 box +mass 1 1.0 + +velocity all create 1.44 87287 loop geom + +pair_style lj/cut 2.5 +pair_coeff 1 1 1.0 1.0 2.5 + +neighbor 0.3 bin +neigh_modify delay 0 every 20 check no + +fix 1 all nve + +variable eng equal pe +variable vy atom vy + +run 100 diff --git a/python/examples/in.gui b/python/examples/in.gui new file mode 100644 index 0000000000..9c88cae484 --- /dev/null +++ b/python/examples/in.gui @@ -0,0 +1,24 @@ +# 3d Lennard-Jones melt + +units lj +dimension 2 +atom_style atomic + +lattice sq2 0.8442 +region box block 0 30 0 15 -0.5 0.5 +create_box 1 box +create_atoms 1 box +mass 1 1.0 + +velocity all create 1.44 87287 loop geom + +pair_style lj/cut 2.5 +pair_coeff 1 1 1.0 1.0 2.5 + +neighbor 0.3 bin +neigh_modify delay 0 every 1 check yes + +timestep 0.003 + +fix 1 all nve +fix 3 all enforce2d diff --git a/python/examples/in.lj b/python/examples/in.lj new file mode 100644 index 0000000000..e69de29bb2 diff --git a/python/examples/in.plot b/python/examples/in.plot new file mode 100644 index 0000000000..9c88cae484 --- /dev/null +++ b/python/examples/in.plot @@ -0,0 +1,24 @@ +# 3d Lennard-Jones melt + +units lj +dimension 2 +atom_style atomic + +lattice sq2 0.8442 +region box block 0 30 0 15 -0.5 0.5 +create_box 1 box +create_atoms 1 box +mass 1 1.0 + +velocity all create 1.44 87287 loop geom + +pair_style lj/cut 2.5 +pair_coeff 1 1 1.0 1.0 2.5 + +neighbor 0.3 bin +neigh_modify delay 0 every 1 check yes + +timestep 0.003 + +fix 1 all nve +fix 3 all enforce2d diff --git a/python/examples/in.simple b/python/examples/in.simple new file mode 100644 index 0000000000..5f6ebdc4f2 --- /dev/null +++ b/python/examples/in.simple @@ -0,0 +1,23 @@ +# 3d Lennard-Jones melt + +units lj +atom_style atomic +atom_modify map array + +lattice fcc 0.8442 +region box block 0 4 0 4 0 4 +create_box 1 box +create_atoms 1 box +mass 1 1.0 + +velocity all create 1.44 87287 loop geom + +pair_style lj/cut 2.5 +pair_coeff 1 1 1.0 1.0 2.5 + +neighbor 0.3 bin +neigh_modify delay 0 every 20 check no + +fix 1 all nve + +run 10 diff --git a/python/examples/in.trivial b/python/examples/in.trivial new file mode 100644 index 0000000000..fa9739d65c --- /dev/null +++ b/python/examples/in.trivial @@ -0,0 +1,22 @@ +# 3d Lennard-Jones melt + +units lj +atom_style atomic + +lattice fcc 0.8442 +region box block 0 10 0 10 0 10 +create_box 1 box +create_atoms 1 box +mass 1 1.0 + +velocity all create 1.44 87287 loop geom + +pair_style lj/cut 2.5 +pair_coeff 1 1 1.0 1.0 2.5 + +neighbor 0.3 bin +neigh_modify delay 0 every 20 check no + +fix 1 all nve + +run 100 diff --git a/python/examples/in.viz b/python/examples/in.viz new file mode 100644 index 0000000000..9c88cae484 --- /dev/null +++ b/python/examples/in.viz @@ -0,0 +1,24 @@ +# 3d Lennard-Jones melt + +units lj +dimension 2 +atom_style atomic + +lattice sq2 0.8442 +region box block 0 30 0 15 -0.5 0.5 +create_box 1 box +create_atoms 1 box +mass 1 1.0 + +velocity all create 1.44 87287 loop geom + +pair_style lj/cut 2.5 +pair_coeff 1 1 1.0 1.0 2.5 + +neighbor 0.3 bin +neigh_modify delay 0 every 1 check yes + +timestep 0.003 + +fix 1 all nve +fix 3 all enforce2d diff --git a/python/examples/log.lammps b/python/examples/log.lammps new file mode 100644 index 0000000000..c8a30e5f53 --- /dev/null +++ b/python/examples/log.lammps @@ -0,0 +1,906 @@ +LAMMPS (31 Oct 2010) +# 3d Lennard-Jones melt + +units lj +dimension 2 +atom_style atomic + +lattice sq2 0.8442 +Lattice spacing in x,y,z = 1.53919 1.53919 1.53919 +region box block 0 30 0 15 -0.5 0.5 +create_box 1 box +Created orthogonal box = (0 0 -0.769595) to (46.1757 23.0879 0.769595) + 1 by 1 by 1 processor grid +create_atoms 1 box +Created 900 atoms +mass 1 1.0 + +velocity all create 1.44 87287 loop geom + +pair_style lj/cut 2.5 +pair_coeff 1 1 1.0 1.0 2.5 + +neighbor 0.3 bin +neigh_modify delay 0 every 1 check yes + +timestep 0.003 + +fix 1 all nve +fix 3 all enforce2d +thermo 100 +fix 2 all langevin 1 1 0.1 12345 +run 100 pre yes post no +Memory usage per processor = 1.67412 Mbytes +Step Temp E_pair E_mol TotEng Press + 0 1.44 -2.6248859 0 -1.1864859 2.0176244 + 100 0.99023288 -2.0005395 0 -1.0114069 5.9914427 +Loop time of 0.04617 on 1 procs for 100 steps with 900 atoms +run 100 pre no post no +Step Temp E_pair E_mol TotEng Press + 100 0.99023288 -2.0005395 0 -1.0114069 5.9914427 + 200 1.0641026 -2.126693 0 -1.0637727 5.3834133 +Loop time of 0.045666 on 1 procs for 100 steps with 900 atoms +run 100 pre no post no +Step Temp E_pair E_mol TotEng Press + 200 1.0641026 -2.126693 0 -1.0637727 5.3834133 + 300 1.0404544 -2.1182885 0 -1.0789901 5.4751694 +Loop time of 0.0448651 on 1 procs for 100 steps with 900 atoms +run 100 pre no post no +Step Temp E_pair E_mol TotEng Press + 300 1.0404544 -2.1182885 0 -1.0789901 5.4751694 + 400 0.9940782 -2.0913772 0 -1.0984035 5.7006406 +Loop time of 0.044888 on 1 procs for 100 steps with 900 atoms +run 100 pre no post no +Step Temp E_pair E_mol TotEng Press + 400 0.9940782 -2.0913772 0 -1.0984035 5.7006406 + 500 1.0533477 -2.1766862 0 -1.1245089 5.1818937 +Loop time of 0.0443962 on 1 procs for 100 steps with 900 atoms +run 100 pre no post no +Step Temp E_pair E_mol TotEng Press + 500 1.0533477 -2.1766862 0 -1.1245089 5.1818937 + 600 1.0371191 -2.1708671 0 -1.1349004 5.2128526 +Loop time of 0.0447521 on 1 procs for 100 steps with 900 atoms +run 100 pre no post no +Step Temp E_pair E_mol TotEng Press + 600 1.0371191 -2.1708671 0 -1.1349004 5.2128526 + 700 1.0178469 -2.1795112 0 -1.1627952 5.0714199 +Loop time of 0.0448291 on 1 procs for 100 steps with 900 atoms +run 100 pre no post no +Step Temp E_pair E_mol TotEng Press + 700 1.0178469 -2.1795112 0 -1.1627952 5.0714199 + 800 1.0220115 -2.1753909 0 -1.1545149 5.1482844 +Loop time of 0.044615 on 1 procs for 100 steps with 900 atoms +run 100 pre no post no +Step Temp E_pair E_mol TotEng Press + 800 1.0220115 -2.1753909 0 -1.1545149 5.1482844 + 900 1.0029029 -2.1674923 0 -1.1657037 5.1338534 +Loop time of 0.0446992 on 1 procs for 100 steps with 900 atoms +run 100 pre no post no +Step Temp E_pair E_mol TotEng Press + 900 1.0029029 -2.1674923 0 -1.1657037 5.1338534 + 1000 1.0400785 -2.1739867 0 -1.1350638 5.1279779 +Loop time of 0.0448501 on 1 procs for 100 steps with 900 atoms +run 100 pre no post no +Step Temp E_pair E_mol TotEng Press + 1000 1.0400785 -2.1739867 0 -1.1350638 5.1279779 + 1100 1.0248051 -2.1763751 0 -1.1527087 5.125636 +Loop time of 0.0447941 on 1 procs for 100 steps with 900 atoms +run 100 pre no post no +Step Temp E_pair E_mol TotEng Press + 1100 1.0248051 -2.1763751 0 -1.1527087 5.125636 + 1200 1.0129247 -2.1960205 0 -1.1842213 4.9280331 +Loop time of 0.0442801 on 1 procs for 100 steps with 900 atoms +run 100 pre no post no +Step Temp E_pair E_mol TotEng Press + 1200 1.0129247 -2.1960205 0 -1.1842213 4.9280331 + 1300 1.0445891 -2.138322 0 -1.0948936 5.358069 +Loop time of 0.044374 on 1 procs for 100 steps with 900 atoms +run 100 pre no post no +Step Temp E_pair E_mol TotEng Press + 1300 1.0445891 -2.138322 0 -1.0948936 5.358069 + 1400 1.0262085 -2.1210161 0 -1.0959478 5.5887816 +Loop time of 0.044693 on 1 procs for 100 steps with 900 atoms +run 100 pre no post no +Step Temp E_pair E_mol TotEng Press + 1400 1.0262085 -2.1210161 0 -1.0959478 5.5887816 + 1500 0.9862776 -2.1287449 0 -1.1435631 5.4408715 +Loop time of 0.0446811 on 1 procs for 100 steps with 900 atoms +run 100 pre no post no +Step Temp E_pair E_mol TotEng Press + 1500 0.9862776 -2.1287449 0 -1.1435631 5.4408715 + 1600 1.0008092 -2.1637199 0 -1.1640227 5.2714443 +Loop time of 0.0442421 on 1 procs for 100 steps with 900 atoms +run 100 pre no post no +Step Temp E_pair E_mol TotEng Press + 1600 1.0008092 -2.1637199 0 -1.1640227 5.2714443 + 1700 1.0258515 -2.1288392 0 -1.1041276 5.4565262 +Loop time of 0.0446551 on 1 procs for 100 steps with 900 atoms +run 100 pre no post no +Step Temp E_pair E_mol TotEng Press + 1700 1.0258515 -2.1288392 0 -1.1041276 5.4565262 + 1800 0.99606403 -2.1881474 0 -1.1931901 5.0145909 +Loop time of 0.0446119 on 1 procs for 100 steps with 900 atoms +run 100 pre no post no +Step Temp E_pair E_mol TotEng Press + 1800 0.99606403 -2.1881474 0 -1.1931901 5.0145909 + 1900 0.98665733 -2.1993707 0 -1.2138096 4.9672864 +Loop time of 0.044836 on 1 procs for 100 steps with 900 atoms +run 100 pre no post no +Step Temp E_pair E_mol TotEng Press + 1900 0.98665733 -2.1993707 0 -1.2138096 4.9672864 + 2000 1.0155716 -2.1239518 0 -1.1095086 5.4592445 +Loop time of 0.044791 on 1 procs for 100 steps with 900 atoms +run 100 pre no post no +Step Temp E_pair E_mol TotEng Press + 2000 1.0155716 -2.1239518 0 -1.1095086 5.4592445 + 2100 0.9688746 -2.1712503 0 -1.2034522 5.1417376 +Loop time of 0.0446949 on 1 procs for 100 steps with 900 atoms +run 100 pre no post no +Step Temp E_pair E_mol TotEng Press + 2100 0.9688746 -2.1712503 0 -1.2034522 5.1417376 + 2200 1.0162371 -2.1732231 0 -1.1581152 5.1108471 +Loop time of 0.0448599 on 1 procs for 100 steps with 900 atoms +run 100 pre no post no +Step Temp E_pair E_mol TotEng Press + 2200 1.0162371 -2.1732231 0 -1.1581152 5.1108471 + 2300 1.0245646 -2.1552369 0 -1.1318107 5.3316927 +Loop time of 0.0446439 on 1 procs for 100 steps with 900 atoms +run 100 pre no post no +Step Temp E_pair E_mol TotEng Press + 2300 1.0245646 -2.1552369 0 -1.1318107 5.3316927 + 2400 1.0532627 -2.168732 0 -1.1166396 5.2318175 +Loop time of 0.044234 on 1 procs for 100 steps with 900 atoms +run 100 pre no post no +Step Temp E_pair E_mol TotEng Press + 2400 1.0532627 -2.168732 0 -1.1166396 5.2318175 + 2500 1.0110981 -2.1283041 0 -1.1183294 5.470307 +Loop time of 0.044275 on 1 procs for 100 steps with 900 atoms +run 100 pre no post no +Step Temp E_pair E_mol TotEng Press + 2500 1.0110981 -2.1283041 0 -1.1183294 5.470307 + 2600 0.99966902 -2.1728149 0 -1.1742566 5.1412453 +Loop time of 0.044683 on 1 procs for 100 steps with 900 atoms +run 100 pre no post no +Step Temp E_pair E_mol TotEng Press + 2600 0.99966902 -2.1728149 0 -1.1742566 5.1412453 + 2700 0.99239191 -2.1637551 0 -1.1724659 5.2165707 +Loop time of 0.0447261 on 1 procs for 100 steps with 900 atoms +run 100 pre no post no +Step Temp E_pair E_mol TotEng Press + 2700 0.99239191 -2.1637551 0 -1.1724659 5.2165707 + 2800 1.0301978 -2.1548308 0 -1.1257777 5.2638595 +Loop time of 0.044596 on 1 procs for 100 steps with 900 atoms +fix 2 all langevin 5 5 0.1 12346 +run 100 pre yes post no +Memory usage per processor = 1.67412 Mbytes +Step Temp E_pair E_mol TotEng Press + 2800 1.0227564 -2.1548308 0 -1.1332108 5.2575845 + 2900 5.0485455 -0.28277455 0 4.7601615 20.746781 +Loop time of 0.0470662 on 1 procs for 100 steps with 900 atoms +run 100 pre no post no +Step Temp E_pair E_mol TotEng Press + 2900 5.0485455 -0.28277455 0 4.7601615 20.746781 + 3000 4.8862233 -0.017789501 0 4.8630047 22.592756 +Loop time of 0.048033 on 1 procs for 100 steps with 900 atoms +run 100 pre no post no +Step Temp E_pair E_mol TotEng Press + 3000 4.8862233 -0.017789501 0 4.8630047 22.592756 + 3100 4.9044567 -0.47959211 0 4.4194152 19.583353 +Loop time of 0.0480561 on 1 procs for 100 steps with 900 atoms +run 100 pre no post no +Step Temp E_pair E_mol TotEng Press + 3100 4.9044567 -0.47959211 0 4.4194152 19.583353 + 3200 5.0879915 -0.14346616 0 4.938872 21.75943 +Loop time of 0.047853 on 1 procs for 100 steps with 900 atoms +run 100 pre no post no +Step Temp E_pair E_mol TotEng Press + 3200 5.0879915 -0.14346616 0 4.938872 21.75943 + 3300 4.8191113 -0.12521782 0 4.688539 21.66605 +Loop time of 0.0484211 on 1 procs for 100 steps with 900 atoms +run 100 pre no post no +Step Temp E_pair E_mol TotEng Press + 3300 4.8191113 -0.12521782 0 4.688539 21.66605 + 3400 5.1800461 -0.14238575 0 5.0319048 21.79606 +Loop time of 0.0481231 on 1 procs for 100 steps with 900 atoms +run 100 pre no post no +Step Temp E_pair E_mol TotEng Press + 3400 5.1800461 -0.14238575 0 5.0319048 21.79606 + 3500 5.0607641 -0.16081652 0 4.8943245 21.486658 +Loop time of 0.0485091 on 1 procs for 100 steps with 900 atoms +run 100 pre no post no +Step Temp E_pair E_mol TotEng Press + 3500 5.0607641 -0.16081652 0 4.8943245 21.486658 + 3600 5.3373422 -0.41828575 0 4.9131261 20.14676 +Loop time of 0.0482152 on 1 procs for 100 steps with 900 atoms +run 100 pre no post no +Step Temp E_pair E_mol TotEng Press + 3600 5.3373422 -0.41828575 0 4.9131261 20.14676 + 3700 5.3202661 -0.30023033 0 5.0141243 20.855754 +Loop time of 0.0481288 on 1 procs for 100 steps with 900 atoms +run 100 pre no post no +Step Temp E_pair E_mol TotEng Press + 3700 5.3202661 -0.30023033 0 5.0141243 20.855754 + 3800 4.8588585 -0.10399834 0 4.7494614 21.771522 +Loop time of 0.0481648 on 1 procs for 100 steps with 900 atoms +run 100 pre no post no +Step Temp E_pair E_mol TotEng Press + 3800 4.8588585 -0.10399834 0 4.7494614 21.771522 + 3900 4.9193657 0.060661037 0 4.9745608 22.688265 +Loop time of 0.0485899 on 1 procs for 100 steps with 900 atoms +run 100 pre no post no +Step Temp E_pair E_mol TotEng Press + 3900 4.9193657 0.060661037 0 4.9745608 22.688265 + 4000 4.9689317 -0.135649 0 4.8277616 21.630946 +Loop time of 0.0486522 on 1 procs for 100 steps with 900 atoms +run 100 pre no post no +Step Temp E_pair E_mol TotEng Press + 4000 4.9689317 -0.135649 0 4.8277616 21.630946 + 4100 5.0796995 -0.42489129 0 4.6491641 19.912881 +Loop time of 0.047446 on 1 procs for 100 steps with 900 atoms +run 100 pre no post no +Step Temp E_pair E_mol TotEng Press + 4100 5.0796995 -0.42489129 0 4.6491641 19.912881 + 4200 4.9121803 -0.26921478 0 4.6375075 20.682997 +Loop time of 0.047961 on 1 procs for 100 steps with 900 atoms +run 100 pre no post no +Step Temp E_pair E_mol TotEng Press + 4200 4.9121803 -0.26921478 0 4.6375075 20.682997 + 4300 5.2480878 -0.25203366 0 4.9902229 21.055825 +Loop time of 0.048105 on 1 procs for 100 steps with 900 atoms +run 100 pre no post no +Step Temp E_pair E_mol TotEng Press + 4300 5.2480878 -0.25203366 0 4.9902229 21.055825 + 4400 5.1754321 -0.1379211 0 5.0317606 21.721423 +Loop time of 0.0476389 on 1 procs for 100 steps with 900 atoms +run 100 pre no post no +Step Temp E_pair E_mol TotEng Press + 4400 5.1754321 -0.1379211 0 5.0317606 21.721423 + 4500 4.8163772 -0.13665697 0 4.6743687 21.338168 +Loop time of 0.048456 on 1 procs for 100 steps with 900 atoms +run 100 pre no post no +Step Temp E_pair E_mol TotEng Press + 4500 4.8163772 -0.13665697 0 4.6743687 21.338168 + 4600 5.1321419 -0.30134093 0 4.8250986 20.640859 +Loop time of 0.048126 on 1 procs for 100 steps with 900 atoms +run 100 pre no post no +Step Temp E_pair E_mol TotEng Press + 4600 5.1321419 -0.30134093 0 4.8250986 20.640859 + 4700 5.054895 -0.11869342 0 4.930585 21.525844 +Loop time of 0.0480442 on 1 procs for 100 steps with 900 atoms +run 100 pre no post no +Step Temp E_pair E_mol TotEng Press + 4700 5.054895 -0.11869342 0 4.930585 21.525844 + 4800 4.7735419 -0.28788546 0 4.4803525 20.241875 +Loop time of 0.048465 on 1 procs for 100 steps with 900 atoms +run 100 pre no post no +Step Temp E_pair E_mol TotEng Press + 4800 4.7735419 -0.28788546 0 4.4803525 20.241875 + 4900 5.1590733 -0.36928517 0 4.7840559 19.906086 +Loop time of 0.0482421 on 1 procs for 100 steps with 900 atoms +run 100 pre no post no +Step Temp E_pair E_mol TotEng Press + 4900 5.1590733 -0.36928517 0 4.7840559 19.906086 + 5000 5.172434 -0.21783749 0 4.9488493 21.017977 +Loop time of 0.0482352 on 1 procs for 100 steps with 900 atoms +run 100 pre no post no +Step Temp E_pair E_mol TotEng Press + 5000 5.172434 -0.21783749 0 4.9488493 21.017977 + 5100 5.3222066 -0.28431917 0 5.0319738 20.654341 +Loop time of 0.0479338 on 1 procs for 100 steps with 900 atoms +run 100 pre no post no +Step Temp E_pair E_mol TotEng Press + 5100 5.3222066 -0.28431917 0 5.0319738 20.654341 + 5200 5.2269814 -0.24120657 0 4.9799671 20.863157 +Loop time of 0.0479622 on 1 procs for 100 steps with 900 atoms +run 100 pre no post no +Step Temp E_pair E_mol TotEng Press + 5200 5.2269814 -0.24120657 0 4.9799671 20.863157 + 5300 4.9617285 -0.25596507 0 4.7002504 20.508874 +Loop time of 0.047616 on 1 procs for 100 steps with 900 atoms +run 100 pre no post no +Step Temp E_pair E_mol TotEng Press + 5300 4.9617285 -0.25596507 0 4.7002504 20.508874 + 5400 5.1034361 -0.23375901 0 4.8640066 20.874041 +Loop time of 0.048399 on 1 procs for 100 steps with 900 atoms +run 100 pre no post no +Step Temp E_pair E_mol TotEng Press + 5400 5.1034361 -0.23375901 0 4.8640066 20.874041 + 5500 5.1033392 -0.023852659 0 5.0738162 21.959054 +Loop time of 0.0481031 on 1 procs for 100 steps with 900 atoms +run 100 pre no post no +Step Temp E_pair E_mol TotEng Press + 5500 5.1033392 -0.023852659 0 5.0738162 21.959054 + 5600 4.9693008 -0.3822577 0 4.5815217 19.745578 +Loop time of 0.0481269 on 1 procs for 100 steps with 900 atoms +run 100 pre no post no +Step Temp E_pair E_mol TotEng Press + 5600 4.9693008 -0.3822577 0 4.5815217 19.745578 + 5700 4.9504059 -0.16971504 0 4.7751904 20.984851 +Loop time of 0.047688 on 1 procs for 100 steps with 900 atoms +run 100 pre no post no +Step Temp E_pair E_mol TotEng Press + 5700 4.9504059 -0.16971504 0 4.7751904 20.984851 + 5800 5.4859652 -0.12581898 0 5.3540507 21.643131 +Loop time of 0.0482688 on 1 procs for 100 steps with 900 atoms +run 100 pre no post no +Step Temp E_pair E_mol TotEng Press + 5800 5.4859652 -0.12581898 0 5.3540507 21.643131 + 5900 5.1041625 -0.24499773 0 4.8534935 20.550919 +Loop time of 0.0485001 on 1 procs for 100 steps with 900 atoms +run 100 pre no post no +Step Temp E_pair E_mol TotEng Press + 5900 5.1041625 -0.24499773 0 4.8534935 20.550919 + 6000 4.7912324 -0.23215106 0 4.5537577 20.385388 +Loop time of 0.0481 on 1 procs for 100 steps with 900 atoms +run 100 pre no post no +Step Temp E_pair E_mol TotEng Press + 6000 4.7912324 -0.23215106 0 4.5537577 20.385388 + 6100 4.8687128 -0.36715522 0 4.4961479 19.616564 +Loop time of 0.047904 on 1 procs for 100 steps with 900 atoms +run 100 pre no post no +Step Temp E_pair E_mol TotEng Press + 6100 4.8687128 -0.36715522 0 4.4961479 19.616564 + 6200 5.1633792 -0.42749917 0 4.7301429 19.411301 +Loop time of 0.0479479 on 1 procs for 100 steps with 900 atoms +run 100 pre no post no +Step Temp E_pair E_mol TotEng Press + 6200 5.1633792 -0.42749917 0 4.7301429 19.411301 + 6300 5.0233319 0.019274316 0 5.0370247 22.144843 +Loop time of 0.048382 on 1 procs for 100 steps with 900 atoms +run 100 pre no post no +Step Temp E_pair E_mol TotEng Press + 6300 5.0233319 0.019274316 0 5.0370247 22.144843 + 6400 4.9283818 -0.080911928 0 4.8419939 21.398618 +Loop time of 0.048492 on 1 procs for 100 steps with 900 atoms +run 100 pre no post no +Step Temp E_pair E_mol TotEng Press + 6400 4.9283818 -0.080911928 0 4.8419939 21.398618 + 6500 5.0852451 -0.29922243 0 4.7803724 20.191962 +Loop time of 0.0481842 on 1 procs for 100 steps with 900 atoms +run 100 pre no post no +Step Temp E_pair E_mol TotEng Press + 6500 5.0852451 -0.29922243 0 4.7803724 20.191962 + 6600 4.9084541 -0.25046783 0 4.6525324 20.294621 +Loop time of 0.048028 on 1 procs for 100 steps with 900 atoms +run 100 pre no post no +Step Temp E_pair E_mol TotEng Press + 6600 4.9084541 -0.25046783 0 4.6525324 20.294621 + 6700 5.158646 -0.24156185 0 4.9113523 20.674148 +Loop time of 0.048152 on 1 procs for 100 steps with 900 atoms +run 100 pre no post no +Step Temp E_pair E_mol TotEng Press + 6700 5.158646 -0.24156185 0 4.9113523 20.674148 + 6800 5.0800581 -0.26396878 0 4.8104448 20.254552 +Loop time of 0.0481989 on 1 procs for 100 steps with 900 atoms +run 100 pre no post no +Step Temp E_pair E_mol TotEng Press + 6800 5.0800581 -0.26396878 0 4.8104448 20.254552 + 6900 5.0775502 -0.39390401 0 4.6780045 19.376358 +Loop time of 0.0482252 on 1 procs for 100 steps with 900 atoms +run 100 pre no post no +Step Temp E_pair E_mol TotEng Press + 6900 5.0775502 -0.39390401 0 4.6780045 19.376358 + 7000 4.9478484 -0.14828675 0 4.7940641 20.564728 +Loop time of 0.0481858 on 1 procs for 100 steps with 900 atoms +run 100 pre no post no +Step Temp E_pair E_mol TotEng Press + 7000 4.9478484 -0.14828675 0 4.7940641 20.564728 + 7100 5.1059459 -0.097059214 0 5.0032134 21.193873 +Loop time of 0.0478609 on 1 procs for 100 steps with 900 atoms +run 100 pre no post no +Step Temp E_pair E_mol TotEng Press + 7100 5.1059459 -0.097059214 0 5.0032134 21.193873 + 7200 5.1752568 -0.37658971 0 4.7929168 19.495134 +Loop time of 0.04796 on 1 procs for 100 steps with 900 atoms +run 100 pre no post no +Step Temp E_pair E_mol TotEng Press + 7200 5.1752568 -0.37658971 0 4.7929168 19.495134 + 7300 5.0979223 -0.059859304 0 5.0323987 21.118636 +Loop time of 0.0480049 on 1 procs for 100 steps with 900 atoms +run 100 pre no post no +Step Temp E_pair E_mol TotEng Press + 7300 5.0979223 -0.059859304 0 5.0323987 21.118636 + 7400 4.980172 -0.27439819 0 4.7002403 19.959942 +Loop time of 0.047915 on 1 procs for 100 steps with 900 atoms +run 100 pre no post no +Step Temp E_pair E_mol TotEng Press + 7400 4.980172 -0.27439819 0 4.7002403 19.959942 + 7500 5.0625014 -0.1884279 0 4.8684485 20.652757 +Loop time of 0.0479548 on 1 procs for 100 steps with 900 atoms +run 100 pre no post no +Step Temp E_pair E_mol TotEng Press + 7500 5.0625014 -0.1884279 0 4.8684485 20.652757 + 7600 4.950022 -0.38785668 0 4.5566653 19.255054 +Loop time of 0.047981 on 1 procs for 100 steps with 900 atoms +run 100 pre no post no +Step Temp E_pair E_mol TotEng Press + 7600 4.950022 -0.38785668 0 4.5566653 19.255054 + 7700 5.038467 0.010501622 0 5.0433703 21.612637 +Loop time of 0.0483868 on 1 procs for 100 steps with 900 atoms +run 100 pre no post no +Step Temp E_pair E_mol TotEng Press + 7700 5.038467 0.010501622 0 5.0433703 21.612637 + 7800 5.0722091 -0.23929391 0 4.8272794 20.163109 +Loop time of 0.0484369 on 1 procs for 100 steps with 900 atoms +run 100 pre no post no +Step Temp E_pair E_mol TotEng Press + 7800 5.0722091 -0.23929391 0 4.8272794 20.163109 + 7900 4.8266996 -0.26738065 0 4.5539559 19.950113 +Loop time of 0.0484252 on 1 procs for 100 steps with 900 atoms +run 100 pre no post no +Step Temp E_pair E_mol TotEng Press + 7900 4.8266996 -0.26738065 0 4.5539559 19.950113 + 8000 5.228568 -0.22189741 0 5.0008611 20.340648 +Loop time of 0.0480399 on 1 procs for 100 steps with 900 atoms +run 100 pre no post no +Step Temp E_pair E_mol TotEng Press + 8000 5.228568 -0.22189741 0 5.0008611 20.340648 + 8100 5.0226139 -0.32338281 0 4.6936504 19.597461 +Loop time of 0.0477419 on 1 procs for 100 steps with 900 atoms +run 100 pre no post no +Step Temp E_pair E_mol TotEng Press + 8100 5.0226139 -0.32338281 0 4.6936504 19.597461 + 8200 4.9067102 -0.1492661 0 4.7519922 20.325845 +Loop time of 0.0479119 on 1 procs for 100 steps with 900 atoms +run 100 pre no post no +Step Temp E_pair E_mol TotEng Press + 8200 4.9067102 -0.1492661 0 4.7519922 20.325845 + 8300 5.1266856 -0.22929386 0 4.8916954 20.201339 +Loop time of 0.0479472 on 1 procs for 100 steps with 900 atoms +run 100 pre no post no +Step Temp E_pair E_mol TotEng Press + 8300 5.1266856 -0.22929386 0 4.8916954 20.201339 + 8400 5.0589881 -0.45921427 0 4.5941528 18.70978 +Loop time of 0.048357 on 1 procs for 100 steps with 900 atoms +run 100 pre no post no +Step Temp E_pair E_mol TotEng Press + 8400 5.0589881 -0.45921427 0 4.5941528 18.70978 + 8500 4.8413482 -0.19493887 0 4.64103 19.954031 +Loop time of 0.0475209 on 1 procs for 100 steps with 900 atoms +run 100 pre no post no +Step Temp E_pair E_mol TotEng Press + 8500 4.8413482 -0.19493887 0 4.64103 19.954031 + 8600 5.0908293 -0.54321499 0 4.5419579 18.109321 +Loop time of 0.0483789 on 1 procs for 100 steps with 900 atoms +run 100 pre no post no +Step Temp E_pair E_mol TotEng Press + 8600 5.0908293 -0.54321499 0 4.5419579 18.109321 + 8700 5.1796752 -0.2916707 0 4.8822493 19.664654 +Loop time of 0.048032 on 1 procs for 100 steps with 900 atoms +run 100 pre no post no +Step Temp E_pair E_mol TotEng Press + 8700 5.1796752 -0.2916707 0 4.8822493 19.664654 + 8800 5.2539084 -0.45997392 0 4.7880968 18.655485 +Loop time of 0.0482469 on 1 procs for 100 steps with 900 atoms +run 100 pre no post no +Step Temp E_pair E_mol TotEng Press + 8800 5.2539084 -0.45997392 0 4.7880968 18.655485 + 8900 5.1024558 -0.23247495 0 4.8643114 19.86178 +Loop time of 0.0481958 on 1 procs for 100 steps with 900 atoms +run 100 pre no post no +Step Temp E_pair E_mol TotEng Press + 8900 5.1024558 -0.23247495 0 4.8643114 19.86178 + 9000 5.1823165 -0.39916589 0 4.7773924 18.971835 +Loop time of 0.0481231 on 1 procs for 100 steps with 900 atoms +run 100 pre no post no +Step Temp E_pair E_mol TotEng Press + 9000 5.1823165 -0.39916589 0 4.7773924 18.971835 + 9100 4.8447578 -0.36539964 0 4.4739751 18.749222 +Loop time of 0.0482869 on 1 procs for 100 steps with 900 atoms +run 100 pre no post no +Step Temp E_pair E_mol TotEng Press + 9100 4.8447578 -0.36539964 0 4.4739751 18.749222 + 9200 4.7762145 -0.41137254 0 4.3595351 18.283329 +Loop time of 0.0483339 on 1 procs for 100 steps with 900 atoms +run 100 pre no post no +Step Temp E_pair E_mol TotEng Press + 9200 4.7762145 -0.41137254 0 4.3595351 18.283329 + 9300 5.0833877 -0.32352098 0 4.7542185 19.215267 +Loop time of 0.0479991 on 1 procs for 100 steps with 900 atoms +run 100 pre no post no +Step Temp E_pair E_mol TotEng Press + 9300 5.0833877 -0.32352098 0 4.7542185 19.215267 + 9400 5.0177659 -0.28646064 0 4.72573 19.284236 +Loop time of 0.048521 on 1 procs for 100 steps with 900 atoms +run 100 pre no post no +Step Temp E_pair E_mol TotEng Press + 9400 5.0177659 -0.28646064 0 4.72573 19.284236 + 9500 5.094124 -0.2545962 0 4.8338677 19.49513 +Loop time of 0.0480311 on 1 procs for 100 steps with 900 atoms +run 100 pre no post no +Step Temp E_pair E_mol TotEng Press + 9500 5.094124 -0.2545962 0 4.8338677 19.49513 + 9600 5.0240503 -0.35703989 0 4.6614281 18.851793 +Loop time of 0.048084 on 1 procs for 100 steps with 900 atoms +run 100 pre no post no +Step Temp E_pair E_mol TotEng Press + 9600 5.0240503 -0.35703989 0 4.6614281 18.851793 + 9700 4.7489492 -0.25367959 0 4.489993 19.113672 +Loop time of 0.047981 on 1 procs for 100 steps with 900 atoms +run 100 pre no post no +Step Temp E_pair E_mol TotEng Press + 9700 4.7489492 -0.25367959 0 4.489993 19.113672 + 9800 4.9909149 -0.34688365 0 4.6384858 18.781061 +Loop time of 0.0480509 on 1 procs for 100 steps with 900 atoms +run 100 pre no post no +Step Temp E_pair E_mol TotEng Press + 9800 4.9909149 -0.34688365 0 4.6384858 18.781061 + 9900 5.2358495 -0.20088636 0 5.0291455 19.968679 +Loop time of 0.0480289 on 1 procs for 100 steps with 900 atoms +run 100 pre no post no +Step Temp E_pair E_mol TotEng Press + 9900 5.2358495 -0.20088636 0 5.0291455 19.968679 + 10000 5.0336433 -0.30070619 0 4.7273442 19.288936 +Loop time of 0.0481491 on 1 procs for 100 steps with 900 atoms +run 100 pre no post no +Step Temp E_pair E_mol TotEng Press + 10000 5.0336433 -0.30070619 0 4.7273442 19.288936 + 10100 4.7368529 -0.43300505 0 4.2985846 18.223831 +Loop time of 0.047874 on 1 procs for 100 steps with 900 atoms +run 100 pre no post no +Step Temp E_pair E_mol TotEng Press + 10100 4.7368529 -0.43300505 0 4.2985846 18.223831 + 10200 5.1921737 -0.39013066 0 4.796274 18.671159 +Loop time of 0.0475681 on 1 procs for 100 steps with 900 atoms +run 100 pre no post no +Step Temp E_pair E_mol TotEng Press + 10200 5.1921737 -0.39013066 0 4.796274 18.671159 + 10300 4.8489929 -0.26376235 0 4.5798428 19.1258 +Loop time of 0.04794 on 1 procs for 100 steps with 900 atoms +run 100 pre no post no +Step Temp E_pair E_mol TotEng Press + 10300 4.8489929 -0.26376235 0 4.5798428 19.1258 + 10400 4.8344238 -0.39563076 0 4.4334214 18.406798 +Loop time of 0.0483451 on 1 procs for 100 steps with 900 atoms +run 100 pre no post no +Step Temp E_pair E_mol TotEng Press + 10400 4.8344238 -0.39563076 0 4.4334214 18.406798 + 10500 4.7947329 -0.30412973 0 4.4852757 18.895224 +Loop time of 0.047621 on 1 procs for 100 steps with 900 atoms +run 100 pre no post no +Step Temp E_pair E_mol TotEng Press + 10500 4.7947329 -0.30412973 0 4.4852757 18.895224 + 10600 5.1411793 -0.14227144 0 4.9931954 19.939056 +Loop time of 0.048378 on 1 procs for 100 steps with 900 atoms +run 100 pre no post no +Step Temp E_pair E_mol TotEng Press + 10600 5.1411793 -0.14227144 0 4.9931954 19.939056 + 10700 4.6126242 -0.35662101 0 4.250878 18.178539 +Loop time of 0.0480251 on 1 procs for 100 steps with 900 atoms +run 100 pre no post no +Step Temp E_pair E_mol TotEng Press + 10700 4.6126242 -0.35662101 0 4.250878 18.178539 + 10800 5.218192 -0.26093796 0 4.9514561 19.419601 +Loop time of 0.0479641 on 1 procs for 100 steps with 900 atoms +run 100 pre no post no +Step Temp E_pair E_mol TotEng Press + 10800 5.218192 -0.26093796 0 4.9514561 19.419601 + 10900 4.9434504 -0.45628696 0 4.4816707 18.016015 +Loop time of 0.048048 on 1 procs for 100 steps with 900 atoms +run 100 pre no post no +Step Temp E_pair E_mol TotEng Press + 10900 4.9434504 -0.45628696 0 4.4816707 18.016015 + 11000 5.2502606 -0.29191121 0 4.9525157 19.228121 +Loop time of 0.0480599 on 1 procs for 100 steps with 900 atoms +run 100 pre no post no +Step Temp E_pair E_mol TotEng Press + 11000 5.2502606 -0.29191121 0 4.9525157 19.228121 + 11100 4.7598343 -0.37412706 0 4.3804185 18.317182 +Loop time of 0.048358 on 1 procs for 100 steps with 900 atoms +run 100 pre no post no +Step Temp E_pair E_mol TotEng Press + 11100 4.7598343 -0.37412706 0 4.3804185 18.317182 + 11200 5.323248 -0.21800357 0 5.0993297 19.810405 +Loop time of 0.04793 on 1 procs for 100 steps with 900 atoms +run 100 pre no post no +Step Temp E_pair E_mol TotEng Press + 11200 5.323248 -0.21800357 0 5.0993297 19.810405 + 11300 4.7219543 -0.50822681 0 4.2084808 17.364301 +Loop time of 0.0486059 on 1 procs for 100 steps with 900 atoms +run 100 pre no post no +Step Temp E_pair E_mol TotEng Press + 11300 4.7219543 -0.50822681 0 4.2084808 17.364301 + 11400 5.2261295 -0.41060675 0 4.8097159 18.184162 +Loop time of 0.0475581 on 1 procs for 100 steps with 900 atoms +run 100 pre no post no +Step Temp E_pair E_mol TotEng Press + 11400 5.2261295 -0.41060675 0 4.8097159 18.184162 + 11500 5.4192661 -0.41967268 0 4.993572 18.44905 +Loop time of 0.048496 on 1 procs for 100 steps with 900 atoms +run 100 pre no post no +Step Temp E_pair E_mol TotEng Press + 11500 5.4192661 -0.41967268 0 4.993572 18.44905 + 11600 4.9327793 -0.29827109 0 4.6290273 18.678345 +Loop time of 0.0480838 on 1 procs for 100 steps with 900 atoms +run 100 pre no post no +Step Temp E_pair E_mol TotEng Press + 11600 4.9327793 -0.29827109 0 4.6290273 18.678345 + 11700 4.9061719 -0.40278149 0 4.4979391 18.015662 +Loop time of 0.0484679 on 1 procs for 100 steps with 900 atoms +run 100 pre no post no +Step Temp E_pair E_mol TotEng Press + 11700 4.9061719 -0.40278149 0 4.4979391 18.015662 + 11800 5.1217904 -0.38967412 0 4.7264254 18.297707 +Loop time of 0.047833 on 1 procs for 100 steps with 900 atoms +run 100 pre no post no +Step Temp E_pair E_mol TotEng Press + 11800 5.1217904 -0.38967412 0 4.7264254 18.297707 + 11900 4.898336 -0.52304599 0 4.3698475 17.482045 +Loop time of 0.0480511 on 1 procs for 100 steps with 900 atoms +run 100 pre no post no +Step Temp E_pair E_mol TotEng Press + 11900 4.898336 -0.52304599 0 4.3698475 17.482045 + 12000 5.0463402 -0.34458395 0 4.6961492 18.259663 +Loop time of 0.048959 on 1 procs for 100 steps with 900 atoms +run 100 pre no post no +Step Temp E_pair E_mol TotEng Press + 12000 5.0463402 -0.34458395 0 4.6961492 18.259663 + 12100 5.2244507 -0.48448431 0 4.7341614 17.921372 +Loop time of 0.0479109 on 1 procs for 100 steps with 900 atoms +run 100 pre no post no +Step Temp E_pair E_mol TotEng Press + 12100 5.2244507 -0.48448431 0 4.7341614 17.921372 + 12200 4.9182331 -0.38972269 0 4.5230457 18.065869 +Loop time of 0.047838 on 1 procs for 100 steps with 900 atoms +run 100 pre no post no +Step Temp E_pair E_mol TotEng Press + 12200 4.9182331 -0.38972269 0 4.5230457 18.065869 + 12300 5.1236535 -0.50804995 0 4.6099106 17.491542 +Loop time of 0.0480042 on 1 procs for 100 steps with 900 atoms +run 100 pre no post no +Step Temp E_pair E_mol TotEng Press + 12300 5.1236535 -0.50804995 0 4.6099106 17.491542 + 12400 4.9892711 -0.35640526 0 4.6273222 18.09428 +Loop time of 0.04794 on 1 procs for 100 steps with 900 atoms +run 100 pre no post no +Step Temp E_pair E_mol TotEng Press + 12400 4.9892711 -0.35640526 0 4.6273222 18.09428 + 12500 4.9708898 -0.47291494 0 4.4924516 17.463457 +Loop time of 0.0483761 on 1 procs for 100 steps with 900 atoms +run 100 pre no post no +Step Temp E_pair E_mol TotEng Press + 12500 4.9708898 -0.47291494 0 4.4924516 17.463457 + 12600 4.9060868 -0.42653866 0 4.4740969 17.831651 +Loop time of 0.047663 on 1 procs for 100 steps with 900 atoms +run 100 pre no post no +Step Temp E_pair E_mol TotEng Press + 12600 4.9060868 -0.42653866 0 4.4740969 17.831651 + 12700 5.0407771 -0.44712314 0 4.5880531 17.595615 +Loop time of 0.0484769 on 1 procs for 100 steps with 900 atoms +run 100 pre no post no +Step Temp E_pair E_mol TotEng Press + 12700 5.0407771 -0.44712314 0 4.5880531 17.595615 + 12800 5.0365185 -0.55979777 0 4.4711246 16.934523 +Loop time of 0.048032 on 1 procs for 100 steps with 900 atoms +run 100 pre no post no +Step Temp E_pair E_mol TotEng Press + 12800 5.0365185 -0.55979777 0 4.4711246 16.934523 + 12900 5.2486399 -0.34366283 0 4.8991453 18.510825 +Loop time of 0.0479882 on 1 procs for 100 steps with 900 atoms +run 100 pre no post no +Step Temp E_pair E_mol TotEng Press + 12900 5.2486399 -0.34366283 0 4.8991453 18.510825 + 13000 5.1541271 -0.40831419 0 4.7400861 17.984905 +Loop time of 0.0485048 on 1 procs for 100 steps with 900 atoms +run 100 pre no post no +Step Temp E_pair E_mol TotEng Press + 13000 5.1541271 -0.40831419 0 4.7400861 17.984905 + 13100 4.9252416 -0.37543615 0 4.544333 18.047342 +Loop time of 0.0478559 on 1 procs for 100 steps with 900 atoms +run 100 pre no post no +Step Temp E_pair E_mol TotEng Press + 13100 4.9252416 -0.37543615 0 4.544333 18.047342 + 13200 4.7618949 -0.57797993 0 4.178624 16.678145 +Loop time of 0.047977 on 1 procs for 100 steps with 900 atoms +run 100 pre no post no +Step Temp E_pair E_mol TotEng Press + 13200 4.7618949 -0.57797993 0 4.178624 16.678145 + 13300 5.0754762 -0.494503 0 4.5753338 17.274009 +Loop time of 0.04761 on 1 procs for 100 steps with 900 atoms +run 100 pre no post no +Step Temp E_pair E_mol TotEng Press + 13300 5.0754762 -0.494503 0 4.5753338 17.274009 + 13400 4.7373562 -0.26050743 0 4.471585 18.244677 +Loop time of 0.048439 on 1 procs for 100 steps with 900 atoms +run 100 pre no post no +Step Temp E_pair E_mol TotEng Press + 13400 4.7373562 -0.26050743 0 4.471585 18.244677 + 13500 5.2850198 -0.52753121 0 4.7516163 17.30467 +Loop time of 0.0479879 on 1 procs for 100 steps with 900 atoms +run 100 pre no post no +Step Temp E_pair E_mol TotEng Press + 13500 5.2850198 -0.52753121 0 4.7516163 17.30467 + 13600 4.962689 -0.43162132 0 4.5255536 17.440072 +Loop time of 0.048393 on 1 procs for 100 steps with 900 atoms +run 100 pre no post no +Step Temp E_pair E_mol TotEng Press + 13600 4.962689 -0.43162132 0 4.5255536 17.440072 + 13700 5.0155844 -0.51839355 0 4.491618 16.999978 +Loop time of 0.048249 on 1 procs for 100 steps with 900 atoms +run 100 pre no post no +Step Temp E_pair E_mol TotEng Press + 13700 5.0155844 -0.51839355 0 4.491618 16.999978 + 13800 5.1030716 -0.54350825 0 4.5538933 17.135096 +Loop time of 0.047961 on 1 procs for 100 steps with 900 atoms +run 100 pre no post no +Step Temp E_pair E_mol TotEng Press + 13800 5.1030716 -0.54350825 0 4.5538933 17.135096 + 13900 4.9484073 -0.32896034 0 4.6139487 17.994498 +Loop time of 0.0479271 on 1 procs for 100 steps with 900 atoms +run 100 pre no post no +Step Temp E_pair E_mol TotEng Press + 13900 4.9484073 -0.32896034 0 4.6139487 17.994498 + 14000 5.1556711 -0.39450883 0 4.7554337 18.024111 +Loop time of 0.0485859 on 1 procs for 100 steps with 900 atoms +run 100 pre no post no +Step Temp E_pair E_mol TotEng Press + 14000 5.1556711 -0.39450883 0 4.7554337 18.024111 + 14100 5.3411027 -0.33237876 0 5.0027894 18.189551 +Loop time of 0.047961 on 1 procs for 100 steps with 900 atoms +run 100 pre no post no +Step Temp E_pair E_mol TotEng Press + 14100 5.3411027 -0.33237876 0 5.0027894 18.189551 + 14200 5.2775385 -0.20079896 0 5.0708756 18.794478 +Loop time of 0.0479829 on 1 procs for 100 steps with 900 atoms +run 100 pre no post no +Step Temp E_pair E_mol TotEng Press + 14200 5.2775385 -0.20079896 0 5.0708756 18.794478 + 14300 4.9698674 -0.54120662 0 4.4231387 16.769686 +Loop time of 0.048353 on 1 procs for 100 steps with 900 atoms +fix 2 all langevin 3 3 0.1 12347 +run 100 pre yes post no +Memory usage per processor = 1.67412 Mbytes +Step Temp E_pair E_mol TotEng Press + 14300 4.9315283 -0.54120662 0 4.3848423 16.737356 + 14400 3.0737982 -1.157576 0 1.9128069 11.390762 +Loop time of 0.0469351 on 1 procs for 100 steps with 900 atoms +run 100 pre no post no +Step Temp E_pair E_mol TotEng Press + 14400 3.0737982 -1.157576 0 1.9128069 11.390762 + 14500 3.0997503 -1.1613569 0 1.9349492 11.746787 +Loop time of 0.046298 on 1 procs for 100 steps with 900 atoms +run 100 pre no post no +Step Temp E_pair E_mol TotEng Press + 14500 3.0997503 -1.1613569 0 1.9349492 11.746787 + 14600 3.0089981 -1.2950599 0 1.7105948 10.694304 +Loop time of 0.0462 on 1 procs for 100 steps with 900 atoms +run 100 pre no post no +Step Temp E_pair E_mol TotEng Press + 14600 3.0089981 -1.2950599 0 1.7105948 10.694304 + 14700 3.0731576 -1.1944523 0 1.8752907 11.177629 +Loop time of 0.0471361 on 1 procs for 100 steps with 900 atoms +run 100 pre no post no +Step Temp E_pair E_mol TotEng Press + 14700 3.0731576 -1.1944523 0 1.8752907 11.177629 + 14800 3.0790691 -1.2454277 0 1.8302202 11.050986 +Loop time of 0.0467501 on 1 procs for 100 steps with 900 atoms +run 100 pre no post no +Step Temp E_pair E_mol TotEng Press + 14800 3.0790691 -1.2454277 0 1.8302202 11.050986 + 14900 2.9875216 -1.2120329 0 1.7721692 11.196195 +Loop time of 0.046731 on 1 procs for 100 steps with 900 atoms +run 100 pre no post no +Step Temp E_pair E_mol TotEng Press + 14900 2.9875216 -1.2120329 0 1.7721692 11.196195 + 15000 3.0023895 -1.1473321 0 1.8517215 11.450509 +Loop time of 0.046834 on 1 procs for 100 steps with 900 atoms +run 100 pre no post no +Step Temp E_pair E_mol TotEng Press + 15000 3.0023895 -1.1473321 0 1.8517215 11.450509 + 15100 3.1095672 -1.186375 0 1.9197372 11.392239 +Loop time of 0.046701 on 1 procs for 100 steps with 900 atoms +run 100 pre no post no +Step Temp E_pair E_mol TotEng Press + 15100 3.1095672 -1.186375 0 1.9197372 11.392239 + 15200 3.0788594 -1.1581741 0 1.9172643 11.434699 +Loop time of 0.0466931 on 1 procs for 100 steps with 900 atoms +run 100 pre no post no +Step Temp E_pair E_mol TotEng Press + 15200 3.0788594 -1.1581741 0 1.9172643 11.434699 + 15300 3.097962 -1.2142726 0 1.8802473 11.060786 +Loop time of 0.0467482 on 1 procs for 100 steps with 900 atoms +run 100 pre no post no +Step Temp E_pair E_mol TotEng Press + 15300 3.097962 -1.2142726 0 1.8802473 11.060786 + 15400 3.0085051 -1.2068705 0 1.7982919 11.125663 +Loop time of 0.046351 on 1 procs for 100 steps with 900 atoms +run 100 pre no post no +Step Temp E_pair E_mol TotEng Press + 15400 3.0085051 -1.2068705 0 1.7982919 11.125663 + 15500 2.91789 -1.2471287 0 1.6675192 10.908068 +Loop time of 0.0467601 on 1 procs for 100 steps with 900 atoms +run 100 pre no post no +Step Temp E_pair E_mol TotEng Press + 15500 2.91789 -1.2471287 0 1.6675192 10.908068 + 15600 3.1592422 -1.320622 0 1.83511 10.315596 +Loop time of 0.0463359 on 1 procs for 100 steps with 900 atoms +run 100 pre no post no +Step Temp E_pair E_mol TotEng Press + 15600 3.1592422 -1.320622 0 1.83511 10.315596 + 15700 3.1361836 -1.16244 0 1.9702589 11.474858 +Loop time of 0.0472879 on 1 procs for 100 steps with 900 atoms +run 100 pre no post no +Step Temp E_pair E_mol TotEng Press + 15700 3.1361836 -1.16244 0 1.9702589 11.474858 + 15800 3.0412387 -1.2445264 0 1.7933331 10.828289 +Loop time of 0.0468159 on 1 procs for 100 steps with 900 atoms +run 100 pre no post no +Step Temp E_pair E_mol TotEng Press + 15800 3.0412387 -1.2445264 0 1.7933331 10.828289 + 15900 3.1318652 -1.1349315 0 1.9934539 11.506448 +Loop time of 0.04685 on 1 procs for 100 steps with 900 atoms +run 100 pre no post no +Step Temp E_pair E_mol TotEng Press + 15900 3.1318652 -1.1349315 0 1.9934539 11.506448 + 16000 3.1012109 -1.2285496 0 1.8692155 11.001584 +Loop time of 0.0467439 on 1 procs for 100 steps with 900 atoms +run 100 pre no post no +Step Temp E_pair E_mol TotEng Press + 16000 3.1012109 -1.2285496 0 1.8692155 11.001584 + 16100 3.0508387 -1.2590515 0 1.7883974 10.774791 +Loop time of 0.046726 on 1 procs for 100 steps with 900 atoms +run 100 pre no post no +Step Temp E_pair E_mol TotEng Press + 16100 3.0508387 -1.2590515 0 1.7883974 10.774791 + 16200 3.1900114 -1.3779639 0 1.8085031 10.051609 +Loop time of 0.047003 on 1 procs for 100 steps with 900 atoms +run 100 pre no post no +Step Temp E_pair E_mol TotEng Press + 16200 3.1900114 -1.3779639 0 1.8085031 10.051609 + 16300 2.9418317 -1.2584583 0 1.6801047 10.537777 +Loop time of 0.0466199 on 1 procs for 100 steps with 900 atoms +run 100 pre no post no +Step Temp E_pair E_mol TotEng Press + 16300 2.9418317 -1.2584583 0 1.6801047 10.537777 + 16400 2.9386319 -1.2427586 0 1.6926082 10.643756 +Loop time of 0.0467391 on 1 procs for 100 steps with 900 atoms +run 100 pre no post no +Step Temp E_pair E_mol TotEng Press + 16400 2.9386319 -1.2427586 0 1.6926082 10.643756 + 16500 2.9129152 -1.0821846 0 1.8274941 11.598824 +Loop time of 0.046731 on 1 procs for 100 steps with 900 atoms +run 100 pre no post no +Step Temp E_pair E_mol TotEng Press + 16500 2.9129152 -1.0821846 0 1.8274941 11.598824 + 16600 3.09857 -1.1510317 0 1.9440954 11.235013 +Loop time of 0.046741 on 1 procs for 100 steps with 900 atoms +run 100 pre no post no +Step Temp E_pair E_mol TotEng Press + 16600 3.09857 -1.1510317 0 1.9440954 11.235013 + 16700 2.8592356 -1.2474647 0 1.6085939 10.595499 +Loop time of 0.046757 on 1 procs for 100 steps with 900 atoms +run 100 pre no post no +Step Temp E_pair E_mol TotEng Press + 16700 2.8592356 -1.2474647 0 1.6085939 10.595499 + 16800 2.930352 -1.2310556 0 1.6960404 10.789207 +Loop time of 0.04637 on 1 procs for 100 steps with 900 atoms +run 100 pre no post no +Step Temp E_pair E_mol TotEng Press + 16800 2.930352 -1.2310556 0 1.6960404 10.789207 + 16900 2.9631052 -1.2091634 0 1.7506495 10.861427 +Loop time of 0.046385 on 1 procs for 100 steps with 900 atoms +run 100 pre no post no +Step Temp E_pair E_mol TotEng Press + 16900 2.9631052 -1.2091634 0 1.7506495 10.861427 + 17000 2.9507403 -1.0775066 0 1.8699552 11.496194 +Loop time of 0.0468409 on 1 procs for 100 steps with 900 atoms +run 100 pre no post no +Step Temp E_pair E_mol TotEng Press + 17000 2.9507403 -1.0775066 0 1.8699552 11.496194 + 17100 3.1429056 -1.1518542 0 1.9875593 11.382163 +Loop time of 0.0466008 on 1 procs for 100 steps with 900 atoms +run 100 pre no post no +Step Temp E_pair E_mol TotEng Press + 17100 3.1429056 -1.1518542 0 1.9875593 11.382163 + 17200 3.0175223 -1.2168934 0 1.7972761 10.888517 +Loop time of 0.0468168 on 1 procs for 100 steps with 900 atoms +run 100 pre no post no +Step Temp E_pair E_mol TotEng Press + 17200 3.0175223 -1.2168934 0 1.7972761 10.888517 + 17300 3.0849889 -1.2641589 0 1.8174022 10.561238 +Loop time of 0.0463462 on 1 procs for 100 steps with 900 atoms +run 100 pre no post no +Step Temp E_pair E_mol TotEng Press + 17300 3.0849889 -1.2641589 0 1.8174022 10.561238 + 17400 2.9622086 -1.2777097 0 1.6812076 10.441674 +Loop time of 0.0466161 on 1 procs for 100 steps with 900 atoms diff --git a/python/examples/pizza/dump.py b/python/examples/pizza/dump.py new file mode 100644 index 0000000000..76e1f3bcee --- /dev/null +++ b/python/examples/pizza/dump.py @@ -0,0 +1,1229 @@ +# 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) + + not needed if dump file is self-describing + +d.tselect.all() select all timesteps +d.tselect.one(N) select only timestep N +d.tselect.none() deselect all timesteps +d.tselect.skip(M) select every Mth step +d.tselect.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.write("file",head,app) write selected steps/atoms to dump file +d.scatter("tmp") write selected steps/atoms to multiple files + + write() can be specified with 2 additional flags + headd = 0/1 for no/yes snapshot header, app = 0/1 for write vs append + 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.setv("type",vector) set a column to a vector of values +d.spread("ke",N,"color") 2nd col = N ints spread over 1st col +d.clone(1000,"color") clone timestep N values to other steps + + 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 + setv() operates on selected timesteps and atoms + if column label does not exist, column is created + values in vector are assigned sequentially to atoms, so may want to sort() + length of vector must match # of selected atoms + 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 +# 12/09, David Hart (SNL): allow use of NumPy or Numeric + +# 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() + +try: + import numpy as np + oldnumeric = False +except: + import Numeric as np + oldnumeric = True + +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 file wasn't self-describing + + if len(self.snaps) == 0: + print "no column assignments made" + elif len(self.names): + print "assigned columns:",self.names2str() + 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:",self.names2str() + elif len(self.snaps[0].atoms[0]) == 8: + self.map(1,"id",2,"type",3,"x",4,"y",5,"z",6,"ix",7,"iy",8,"iz") + print "assigned columns:",self.names2str() + else: + print "no column assignments made" + + # if snapshots are scaled, unscale them + + if (not self.names.has_key("x")) or \ + (not self.names.has_key("y")) or \ + (not self.names.has_key("z")): + print "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],'rb') + 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 + # assign column names if not already done and file is self-describing + # convert xs,xu to x + + def read_snapshot(self,f): + try: + snap = Snap() + item = f.readline() + snap.time = int(f.readline().split()[0]) # just grab 1st field + item = f.readline() + snap.natoms = int(f.readline()) + + snap.aselect = np.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 len(self.names) == 0: + words = item.split()[2:] + if len(words): + for i in range(len(words)): + if words[i] == "xs" or words[i] == "xu": + self.names["x"] = i + elif words[i] == "ys" or words[i] == "yu": + self.names["y"] = i + elif words[i] == "zs" or words[i] == "zu": + self.names["z"] = i + else: self.names[words[i]] = i + + 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) + if oldnumeric: atoms = np.zeros((snap.natoms,ncol),np.Float) + else: atoms = np.zeros((snap.natoms,ncol),np.float) + start = 0 + stop = ncol + for i in xrange(snap.natoms): + atoms[i] = floats[start:stop] + start = stop + stop += ncol + else: atoms = None + snap.atoms = atoms + 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 + + # -------------------------------------------------------------------- + # convert column names assignment to a string, in column order + + def names2str(self): + ncol = len(self.snaps[0].atoms[0]) + pairs = self.names.items() + values = self.names.values() + str = "" + for i in xrange(ncol): + if i in values: str += pairs[values.index(i)][0] + ' ' + return str + + # -------------------------------------------------------------------- + # 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 = np.argsort(ids) + for i in xrange(len(atoms[0])): + atoms[:,i] = np.take(atoms[:,i],ordering) + + # -------------------------------------------------------------------- + # write a single dump file from current selection + + def write(self,file,header=1,append=0): + if len(self.snaps): namestr = self.names2str() + if not append: f = open(file,"w") + else: f = open(file,"a") + for snap in self.snaps: + if not snap.tselect: continue + print snap.time, + sys.stdout.flush() + + if header: + 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",namestr + + 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): + if len(self.snaps): namestr = self.names2str() + 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",namestr + + 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 for all selected snapshots + + 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 + + # -------------------------------------------------------------------- + # set a column value via an input vec for all selected snapshots/atoms + + def setv(self,colname,vec): + print "Setting ..." + if not self.names.has_key(colname): + self.newcolumn(colname) + icol = self.names[colname] + + for snap in self.snaps: + if not snap.tselect: continue + if snap.nselect != len(vec): + raise StandardError,"vec length does not match # of selected atoms" + atoms = snap.atoms + m = 0 + for i in xrange(snap.natoms): + if snap.aselect[i]: + atoms[i][icol] = vec[m] + m += 1 + + # -------------------------------------------------------------------- + # 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 + if oldnumeric: newatoms = np.zeros((snap.natoms,ncol+1),np.Float) + else: newatoms = np.zeros((snap.natoms,ncol+1),np.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/Numpy 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/Numpy 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 + + if oldnumeric: bondlist = np.zeros((nbonds,4),np.Int) + else: bondlist = np.zeros((nbonds,4),np.int) + ints = [abs(int(value)) for value in words] + start = 0 + stop = 4 + for i in 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 --git a/python/examples/pizza/gl.py b/python/examples/pizza/gl.py new file mode 100644 index 0000000000..71db5b63e7 --- /dev/null +++ b/python/examples/pizza/gl.py @@ -0,0 +1,1329 @@ +# 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. + +# gl tool + +oneline = "3d interactive visualization via OpenGL" + +docstr = """ +g = gl(d) create OpenGL display for data in d + + d = atom snapshot object (dump, data) + +g.bg("black") set background color (def = "black") +g.size(N) set image size to NxN +g.size(N,M) set image size to NxM +g.rotate(60,135) view from z theta and azimuthal phi (def = 60,30) +g.shift(x,y) translate by x,y pixels in view window (def = 0,0) +g.zoom(0.5) scale image by factor (def = 1) +g.box(0/1/2) 0/1/2 = none/variable/fixed box +g.box(0/1/2,"green") set box color +g.box(0/1/2,"red",4) set box edge thickness +g.file = "image" file prefix for created images (def = "image") + +g.show(N) show image of snapshot at timestep N + +g.all() make images of all selected snapshots +g.all(P) images of all, start file label at P +g.all(N,M,P) make M images of snapshot N, start label at P + +g.pan(60,135,1.0,40,135,1.5) pan during all() operation +g.pan() no pan during all() (default) + + args = z theta, azimuthal phi, zoom factor at beginning and end + values at each step are interpolated between beginning and end values + +g.select = "$x > %g*3.0" string to pass to d.aselect.test() during all() +g.select = "" no extra aselect (default) + + %g varies from 0.0 to 1.0 from beginning to end of all() + +g.acol(2,"green") set atom colors by atom type (1-N) +g.acol([2,4],["red","blue"]) 1st arg = one type or list of types +g.acol(0,"blue") 2nd arg = one color or list of colors +g.acol(range(20),["red","blue"]) if list lengths unequal, interpolate +g.acol(range(10),"loop") assign colors in loop, randomly ordered + + if 1st arg is 0, set all types to 2nd arg + if list of types has a 0 (e.g. range(10)), +1 is added to each value + interpolate means colors blend smoothly from one value to the next + +g.arad([1,2],[0.5,0.3]) set atom radii, same rules as acol() + +g.bcol() set bond color, same args as acol() +g.brad() set bond thickness, same args as arad() + +g.tcol() set triangle color, same args as acol() +g.tfill() set triangle fill, 0 fill, 1 line, 2 both + +g.lcol() set line color, same args as acol() +g.lrad() set line thickness, same args as arad() + +g.adef() set atom/bond/tri/line properties to default +g.bdef() default = "loop" for colors, 0.45 for radii +g.tdef() default = 0.25 for bond/line thickness +g.ldef() default = 0 fill + + by default 100 types are assigned + if atom/bond/tri/line has type > # defined properties, is an error + +from vizinfo import colors access color list +print colors list defined color names and RGB values +colors["nickname"] = [R,G,B] set new RGB values from 0 to 255 + + 140 pre-defined colors: red, green, blue, purple, yellow, black, white, etc + +Settings specific to gl tool: + +g.q(10) set quality of image (def = 5) +g.axis(0/1) turn xyz axes off/on +g.ortho(0/1) perspective (0) vs orthographic (1) view +g.clip('xlo',0.25) clip in xyz from lo/hi at box fraction (0-1) +g.reload() force all data to be reloaded +g.cache = 0/1 turn off/on GL cache lists (def = on) +theta,phi,x,y,scale,up = g.gview() grab all current view parameters +g.sview(theta,phi,x,y,scale,up) set all view parameters + + data reload is necessary if dump selection is used to change the data + cache lists usually improve graphics performance + gview returns values to use in other commands: + theta,phi are args to rotate() + x,y are args to shift() + scale is arg to zoom() + up is a 3-vector arg to sview() +""" + +# History +# 9/05, Steve Plimpton (SNL): original version + +# ToDo list +# when do aselect with select str while looping N times on same timestep +# would not let you grow # of atoms selected + +# Variables +# ztheta = vertical angle from z-azis of viewpoint +# azphi = azimuthal angle of viewpoint +# xshift,yshift = xy translation of scene (in pixels) +# distance = size of simulation box (largest dim) +# eye = viewpoint distance from center of scene +# file = filename prefix to use for images produced +# boxflag = 0/1/2 for drawing simulation box: none/variable/fixed +# bxcol = color of box +# bxthick = thickness of box lines +# bgcol = color of background +# vizinfo = scene attributes +# center[3] = center point of simulation box +# view[3] = direction towards eye in simulation box (unit vector) +# up[3] = screen up direction in simulation box (unit vector) +# right[3] = screen right direction in simulation box (unit vector) + +# Imports and external programs + +from math import sin,cos,sqrt,pi,acos +from OpenGL.Tk import * +from OpenGL.GLUT import * +import Image +from vizinfo import vizinfo + +# Class definition + +class gl: + +# -------------------------------------------------------------------- + + def __init__(self,data): + self.data = data + self.root = None + self.xpixels = 512 + self.ypixels = 512 + self.ztheta = 60 + self.azphi = 30 + self.scale = 1.0 + self.xshift = self.yshift = 0 + + self.file = "image" + self.boxflag = 0 + self.bxcol = [1,1,0] + self.bxthick = 0.3 + self.bgcol = [0,0,0] + self.labels = [] + self.panflag = 0 + self.select = "" + + self.axisflag = 0 + self.orthoflag = 1 + self.nslices = 5 + self.nstacks = 5 + self.nsides = 10 + self.theta_amplify = 2 + self.shiny = 2 + + self.clipflag = 0 + self.clipxlo = self.clipylo = self.clipzlo = 0.0 + self.clipxhi = self.clipyhi = self.clipzhi = 1.0 + + self.nclist = 0 + self.calllist = [0] # indexed by 1-Ntype, so start with 0 index + self.cache = 1 + self.cachelist = 0 + + self.boxdraw = [] + self.atomdraw = [] + self.bonddraw = [] + self.tridraw = [] + self.linedraw = [] + + self.ready = 0 + self.create_window() + + self.vizinfo = vizinfo() + self.adef() + self.bdef() + self.tdef() + self.ldef() + + self.center = 3*[0] + self.view = 3*[0] + self.up = 3*[0] + self.right = 3*[0] + self.viewupright() + + # -------------------------------------------------------------------- + + def bg(self,color): + from vizinfo import colors + self.bgcol = [colors[color][0]/255.0,colors[color][1]/255.0, + colors[color][2]/255.0] + self.w.tkRedraw() + + # -------------------------------------------------------------------- + + def size(self,xnew,ynew=None): + self.xpixels = xnew + if not ynew: self.ypixels = self.xpixels + else: self.ypixels = ynew + self.create_window() + + # -------------------------------------------------------------------- + + def axis(self,value): + self.axisflag = value + self.cachelist = -self.cachelist + self.w.tkRedraw() + + # -------------------------------------------------------------------- + + def create_window(self): + if self.root: self.root.destroy() + + from __main__ import tkroot + self.root = Toplevel(tkroot) + self.root.title('Pizza.py gl tool') + + self.w = MyOpengl(self.root,width=self.xpixels,height=self.ypixels, + double=1,depth=1) + self.w.pack(expand=YES) +# self.w.pack(expand=YES,fill=BOTH) + + glViewport(0,0,self.xpixels,self.ypixels) + glEnable(GL_LIGHTING); + glEnable(GL_LIGHT0); + glEnable(GL_DEPTH_TEST); + glLightModeli(GL_LIGHT_MODEL_TWO_SIDE,GL_TRUE); + glPolygonMode(GL_FRONT_AND_BACK,GL_FILL) + + self.rtrack = self.xpixels + if self.ypixels > self.xpixels: self.rtrack = self.ypixels + + self.w.redraw = self.redraw + self.w.parent = self + self.w.tkRedraw() + tkroot.update_idletasks() # force window to appear + + # -------------------------------------------------------------------- + + def clip(self,which,value): + if which == "xlo": + self.clipxlo = value + if value > self.clipxhi: self.clipxlo = self.clipxhi + elif which == "xhi": + self.clipxhi = value + if value < self.clipxlo: self.clipxhi = self.clipxlo + elif which == "ylo": + self.clipylo = value + if value > self.clipyhi: self.clipylo = self.clipyhi + elif which == "yhi": + self.clipyhi = value + if value < self.clipylo: self.clipyhi = self.clipylo + elif which == "zlo": + self.clipzlo = value + if value > self.clipzhi: self.clipzlo = self.clipzhi + elif which == "zhi": + self.clipzhi = value + if value < self.clipzlo: self.clipzhi = self.clipzlo + + oldflag = self.clipflag + if self.clipxlo > 0 or self.clipylo > 0 or self.clipzlo > 0 or \ + self.clipxhi < 1 or self.clipyhi < 1 or self.clipzhi < 1: + self.clipflag = 1 + else: self.clipflag = 0 + + if oldflag == 0 and self.clipflag == 0: return + self.cachelist = -self.cachelist + self.w.tkRedraw() + + # -------------------------------------------------------------------- + + def q(self,value): + self.nslices = value + self.nstacks = value + self.make_atom_calllist() + self.cachelist = -self.cachelist + self.w.tkRedraw() + + # -------------------------------------------------------------------- + + def ortho(self,value): + self.orthoflag = value + self.w.tkRedraw() + + # -------------------------------------------------------------------- + # set unit vectors for view,up,right from ztheta,azphi + # assume +z in scene should be up on screen (unless looking down z-axis) + # right = up x view + + def viewupright(self): + self.view[0] = cos(pi*self.azphi/180) * sin(pi*self.ztheta/180) + self.view[1] = sin(pi*self.azphi/180) * sin(pi*self.ztheta/180) + self.view[2] = cos(pi*self.ztheta/180) + + if self.ztheta == 0.0: + self.up[0] = cos(pi*self.azphi/180) + self.up[1] = -sin(pi*self.azphi/180) + self.up[2] = 0.0 + elif self.ztheta == 180.0: + self.up[0] = cos(pi*self.azphi/180) + self.up[1] = sin(pi*self.azphi/180) + self.up[2] = 0.0 + else: + dot = self.view[2] # dot = (0,0,1) . view + self.up[0] = -dot*self.view[0] # up projected onto v = dot * v + self.up[1] = -dot*self.view[1] # up perp to v = up - dot * v + self.up[2] = 1.0 - dot*self.view[2] + + self.up = vecnorm(self.up) + self.right = veccross(self.up,self.view) + + # -------------------------------------------------------------------- + # reset ztheta,azphi and thus view,up.right + # called as function from Pizza.py + + def rotate(self,ztheta,azphi): + self.ztheta = ztheta + self.azphi = azphi + self.viewupright() + self.setview() + self.w.tkRedraw() + + # -------------------------------------------------------------------- + # return all view params to reproduce current display via sview() + + def gview(self): + return self.ztheta,self.azphi,self.xshift,self.yshift,self.scale,self.up + + # -------------------------------------------------------------------- + # set current view, called by user with full set of view params + # up is not settable via any other call, all other params are + + def sview(self,ztheta,azphi,xshift,yshift,scale,up): + self.ztheta = ztheta + self.azphi = azphi + self.xshift = xshift + self.yshift = yshift + self.scale = scale + self.up[0] = up[0] + self.up[1] = up[1] + self.up[2] = up[2] + self.up = vecnorm(self.up) + self.view[0] = cos(pi*self.azphi/180) * sin(pi*self.ztheta/180) + self.view[1] = sin(pi*self.azphi/180) * sin(pi*self.ztheta/180) + self.view[2] = cos(pi*self.ztheta/180) + self.right = veccross(self.up,self.view) + self.setview() + self.w.tkRedraw() + + # -------------------------------------------------------------------- + # rotation triggered by mouse trackball + # project old,new onto unit trackball surf + # rotate view,up around axis of rotation = old x new + # right = up x view + # reset ztheta,azphi from view + + def mouse_rotate(self,xnew,ynew,xold,yold): + + # change y pixels to measure from bottom of window instead of top + + yold = self.ypixels - yold + ynew = self.ypixels - ynew + + # vold = unit vector to (xold,yold) projected onto trackball + # vnew = unit vector to (xnew,ynew) projected onto trackball + # return (no rotation) if either projection point is outside rtrack + + vold = [0,0,0] + vold[0] = xold - (0.5*self.xpixels + self.xshift) + vold[1] = yold - (0.5*self.ypixels + self.yshift) + vold[2] = self.rtrack*self.rtrack - vold[0]*vold[0] - vold[1]*vold[1] + if vold[2] < 0: return + vold[2] = sqrt(vold[2]) + vold = vecnorm(vold) + + vnew = [0,0,0] + vnew[0] = xnew - (0.5*self.xpixels + self.xshift) + vnew[1] = ynew - (0.5*self.ypixels + self.yshift) + vnew[2] = self.rtrack*self.rtrack - vnew[0]*vnew[0] - vnew[1]*vnew[1] + if vnew[2] < 0: return + vnew[2] = sqrt(vnew[2]) + vnew = vecnorm(vnew) + + # rot = trackball rotation axis in screen ref frame = vold x vnew + # theta = angle of rotation = sin(theta) for small theta + # axis = rotation axis in body ref frame described by right,up,view + + rot = veccross(vold,vnew) + theta = sqrt(rot[0]*rot[0] + rot[1]*rot[1] + rot[2]*rot[2]) + theta *= self.theta_amplify + + axis = [0,0,0] + axis[0] = rot[0]*self.right[0] + rot[1]*self.up[0] + rot[2]*self.view[0] + axis[1] = rot[0]*self.right[1] + rot[1]*self.up[1] + rot[2]*self.view[1] + axis[2] = rot[0]*self.right[2] + rot[1]*self.up[2] + rot[2]*self.view[2] + axis = vecnorm(axis) + + # view is changed by (axis x view) scaled by theta + # up is changed by (axis x up) scaled by theta + # force up to be perp to view via up_perp = up - (up . view) view + # right = up x view + + delta = veccross(axis,self.view) + self.view[0] -= theta*delta[0] + self.view[1] -= theta*delta[1] + self.view[2] -= theta*delta[2] + self.view = vecnorm(self.view) + + delta = veccross(axis,self.up) + self.up[0] -= theta*delta[0] + self.up[1] -= theta*delta[1] + self.up[2] -= theta*delta[2] + + dot = vecdot(self.up,self.view) + self.up[0] -= dot*self.view[0] + self.up[1] -= dot*self.view[1] + self.up[2] -= dot*self.view[2] + self.up = vecnorm(self.up) + + self.right = veccross(self.up,self.view) + + # convert new view to ztheta,azphi + + self.ztheta = acos(self.view[2])/pi * 180.0 + if (self.ztheta == 0.0): self.azphi = 0.0 + else: self.azphi = acos(self.view[0]/sin(pi*self.ztheta/180.0))/pi * 180.0 + if self.view[1] < 0: self.azphi = 360.0 - self.azphi + self.setview() + self.w.tkRedraw() + + # -------------------------------------------------------------------- + + def shift(self,x,y): + self.xshift = x; + self.yshift = y; + self.setview() + self.w.tkRedraw() + + # -------------------------------------------------------------------- + + def zoom(self,scale): + self.scale = scale + self.setview() + self.w.tkRedraw() + + # -------------------------------------------------------------------- + # set view params needed by redraw + # input: center = center of box + # distance = size of scene (longest box length) + # scale = zoom factor (1.0 = no zoom) + # xshift,yshift = translation factor in pixels + # view = unit vector from center to viewpoint + # up = unit vector in up direction in scene + # right = unit vector in right direction in scene + # output: eye = distance to view scene from + # xto,yto,zto = point to look to + # xfrom,yfrom,zfrom = point to look from + + def setview(self): + if not self.ready: return # no distance since no scene yet + + self.eye = 3 * self.distance / self.scale + xfactor = 0.5*self.eye*self.xshift/self.xpixels + yfactor = 0.5*self.eye*self.yshift/self.ypixels + + self.xto = self.center[0] - xfactor*self.right[0] - yfactor*self.up[0] + self.yto = self.center[1] - xfactor*self.right[1] - yfactor*self.up[1] + self.zto = self.center[2] - xfactor*self.right[2] - yfactor*self.up[2] + + self.xfrom = self.xto + self.eye*self.view[0] + self.yfrom = self.yto + self.eye*self.view[1] + self.zfrom = self.zto + self.eye*self.view[2] + + # -------------------------------------------------------------------- + # box attributes, also used for triangle lines + + def box(self,*args): + self.boxflag = args[0] + if len(args) > 1: + from vizinfo import colors + self.bxcol = [colors[args[1]][0]/255.0,colors[args[1]][1]/255.0, + colors[args[1]][2]/255.0] + if len(args) > 2: self.bxthick = args[2] + self.cachelist = -self.cachelist + self.w.tkRedraw() + + # -------------------------------------------------------------------- + # grab all selected snapshots from data object + # add GL-specific info to each bond + + def reload(self): + print "Loading data into gl tool ..." + data = self.data + + self.timeframes = [] + self.boxframes = [] + self.atomframes = [] + self.bondframes = [] + self.triframes = [] + self.lineframes = [] + + box = [] + if self.boxflag == 2: box = data.maxbox() + + flag = 0 + while 1: + which,time,flag = data.iterator(flag) + if flag == -1: break + time,boxone,atoms,bonds,tris,lines = data.viz(which) + if self.boxflag < 2: box = boxone + if bonds: self.bonds_augment(bonds) + + self.timeframes.append(time) + self.boxframes.append(box) + self.atomframes.append(atoms) + self.bondframes.append(bonds) + self.triframes.append(tris) + self.lineframes.append(lines) + + print time, + sys.stdout.flush() + print + + self.nframes = len(self.timeframes) + self.distance = compute_distance(self.boxframes[0]) + self.center = compute_center(self.boxframes[0]) + self.ready = 1 + self.setview() + + # -------------------------------------------------------------------- + + def nolabel(self): + self.cachelist = -self.cachelist + self.labels = [] + + # -------------------------------------------------------------------- + # show a single snapshot + # distance from snapshot box or max box for all selected steps + + def show(self,ntime): + data = self.data + which = data.findtime(ntime) + time,box,atoms,bonds,tris,lines = data.viz(which) + if self.boxflag == 2: box = data.maxbox() + self.distance = compute_distance(box) + self.center = compute_center(box) + + if bonds: self.bonds_augment(bonds) + + self.boxdraw = box + self.atomdraw = atoms + self.bonddraw = bonds + self.tridraw = tris + self.linedraw = lines + + self.ready = 1 + self.setview() + self.cachelist = -self.cachelist + self.w.tkRedraw() + self.save() + + # -------------------------------------------------------------------- + + def pan(self,*list): + if len(list) == 0: self.panflag = 0 + else: + self.panflag = 1 + self.ztheta_start = list[0] + self.azphi_start = list[1] + self.scale_start = list[2] + self.ztheta_stop = list[3] + self.azphi_stop = list[4] + self.scale_stop = list[5] + + # -------------------------------------------------------------------- + + def all(self,*list): + data = self.data + if len(list) == 0: + nstart = 0 + ncount = data.nselect + elif len(list) == 1: + nstart = list[0] + ncount = data.nselect + else: + ntime = list[0] + nstart = list[2] + ncount = list[1] + + if self.boxflag == 2: box = data.maxbox() + + # loop over all selected steps + # distance from 1st snapshot box or max box for all selected steps + # recompute box center on 1st step or if panning + + if len(list) <= 1: + + n = nstart + i = flag = 0 + while 1: + which,time,flag = data.iterator(flag) + if flag == -1: break + + fraction = float(i) / (ncount-1) + + if self.select != "": + newstr = self.select % fraction + data.aselect.test(newstr,time) + time,boxone,atoms,bonds,tris,lines = data.viz(which) + + if self.boxflag < 2: box = boxone + if n == nstart: self.distance = compute_distance(box) + + if n < 10: file = self.file + "000" + str(n) + elif n < 100: file = self.file + "00" + str(n) + elif n < 1000: file = self.file + "0" + str(n) + else: file = self.file + str(n) + + if self.panflag: + self.ztheta = self.ztheta_start + \ + fraction*(self.ztheta_stop - self.ztheta_start) + self.azphi = self.azphi_start + \ + fraction*(self.azphi_stop - self.azphi_start) + self.scale = self.scale_start + \ + fraction*(self.scale_stop - self.scale_start) + self.viewupright() + + if n == nstart or self.panflag: self.center = compute_center(box) + + if bonds: self.bonds_augment(bonds) + + self.boxdraw = box + self.atomdraw = atoms + self.bonddraw = bonds + self.tridraw = tris + self.linedraw = lines + + self.ready = 1 + self.setview() + self.cachelist = -self.cachelist + self.w.tkRedraw() + self.save(file) + + print time, + sys.stdout.flush() + i += 1 + n += 1 + + # loop ncount times on same step + # distance from 1st snapshot box or max box for all selected steps + # recompute box center on 1st step or if panning + + else: + + which = data.findtime(ntime) + + n = nstart + for i in range(ncount): + fraction = float(i) / (ncount-1) + + if self.select != "": + newstr = self.select % fraction + data.aselect.test(newstr,ntime) + time,boxone,atoms,bonds,tris,lines = data.viz(which) + + if self.boxflag < 2: box = boxone + if n == nstart: self.distance = compute_distance(box) + + if n < 10: file = self.file + "000" + str(n) + elif n < 100: file = self.file + "00" + str(n) + elif n < 1000: file = self.file + "0" + str(n) + else: file = self.file + str(n) + + if self.panflag: + self.ztheta = self.ztheta_start + \ + fraction*(self.ztheta_stop - self.ztheta_start) + self.azphi = self.azphi_start + \ + fraction*(self.azphi_stop - self.azphi_start) + self.scale = self.scale_start + \ + fraction*(self.scale_stop - self.scale_start) + self.viewupright() + + if n == nstart or self.panflag: self.center = compute_center(box) + + if bonds: self.bonds_augment(bonds) + + self.boxdraw = box + self.atomdraw = atoms + self.bonddraw = bonds + self.tridraw = tris + self.linedraw = lines + + self.ready = 1 + self.setview() + self.cachelist = -self.cachelist + self.w.tkRedraw() + self.save(file) + + print n, + sys.stdout.flush() + n += 1 + + print "\n%d images" % ncount + + # -------------------------------------------------------------------- + + def display(self,index): + self.boxdraw = self.boxframes[index] + self.atomdraw = self.atomframes[index] + self.bonddraw = self.bondframes[index] + self.tridraw = self.triframes[index] + self.linedraw = self.lineframes[index] + + self.ready = 1 + self.cachelist = -self.cachelist + self.w.tkRedraw() + return (self.timeframes[index],len(self.atomdraw)) + + # -------------------------------------------------------------------- + # draw the GL scene + + def redraw(self,o): + # clear window to background color + + glClearColor(self.bgcol[0],self.bgcol[1],self.bgcol[2],0) + glClear(GL_COLOR_BUFFER_BIT | GL_DEPTH_BUFFER_BIT) + + # not ready if no scene yet + + if not self.ready: return + + # set view from eye, distance, 3 lookat vectors (from,to,up) + + glMatrixMode(GL_PROJECTION) + glLoadIdentity() + if self.orthoflag: + glOrtho(-0.25*self.eye,0.25*self.eye,-0.25*self.eye,0.25*self.eye, + self.eye-2*self.distance,self.eye+2*self.distance) + else: + gluPerspective(30.0,1.0,0.01,10000.0) + + glMatrixMode(GL_MODELVIEW) + glLoadIdentity() + gluLookAt(self.xfrom,self.yfrom,self.zfrom,self.xto,self.yto,self.zto, + self.up[0],self.up[1],self.up[2]) + + # draw scene from display list if caching allowed and list hasn't changed + # else redraw and store as new display list if caching allowed + + if self.cache and self.cachelist > 0: glCallList(self.cachelist); + else: + if self.cache: + if self.cachelist < 0: glDeleteLists(-self.cachelist,1) + self.cachelist = glGenLists(1) + glNewList(self.cachelist,GL_COMPILE_AND_EXECUTE) + + # draw box, clip-box, xyz axes, lines + + glDisable(GL_LIGHTING) + + if self.boxflag: + self.draw_box(0) + if self.clipflag: self.draw_box(1) + if self.axisflag: self.draw_axes() + + ncolor = self.vizinfo.nlcolor + for line in self.linedraw: + itype = int(line[1]) + if itype > ncolor: raise StandardError,"line type too big" + red,green,blue = self.vizinfo.lcolor[itype] + glColor3f(red,green,blue) + thick = self.vizinfo.lrad[itype] + glLineWidth(thick) + glBegin(GL_LINES) + glVertex3f(line[2],line[3],line[4]) + glVertex3f(line[5],line[6],line[7]) + glEnd() + + glEnable(GL_LIGHTING) + + # draw non-clipped scene = atoms, bonds, triangles + +# draw atoms as collection of points +# cannot put PointSize inside glBegin +# so probably need to group atoms by type for best performance +# or just allow one radius +# need to scale radius appropriately with box size +# or could leave it at absolute value +# use POINT_SMOOTH to enable anti-aliasing and round points +# multiple timesteps via vcr::play() is still not fast +# caching makes it fast for single frame, but multiple frames is slow +# need to enable clipping + +# if not self.clipflag: +# glDisable(GL_LIGHTING) +# glEnable(GL_POINT_SMOOTH) +# glPointSize(self.vizinfo.arad[int(self.atomdraw[0][1])]) +# glBegin(GL_POINTS) +# for atom in self.atomdraw: +# red,green,blue = self.vizinfo.acolor[int(atom[1])] +# glColor(red,green,blue) +# glVertex3d(atom[2],atom[3],atom[4]) +# glEnd() +# glEnable(GL_LIGHTING) + + if not self.clipflag: + for atom in self.atomdraw: + glTranslatef(atom[2],atom[3],atom[4]); + glCallList(self.calllist[int(atom[1])]); + glTranslatef(-atom[2],-atom[3],-atom[4]); + + if self.bonddraw: + bound = 0.25 * self.distance + ncolor = self.vizinfo.nbcolor + for bond in self.bonddraw: + if bond[10] > bound: continue + itype = int(bond[1]) + if itype > ncolor: raise StandardError,"bond type too big" + red,green,blue = self.vizinfo.bcolor[itype] + rad = self.vizinfo.brad[itype] + glPushMatrix() + glTranslatef(bond[2],bond[3],bond[4]) + glRotatef(bond[11],bond[12],bond[13],0.0) + glMaterialfv(GL_FRONT_AND_BACK,GL_EMISSION,[red,green,blue,1.0]); + glMaterialf(GL_FRONT_AND_BACK,GL_SHININESS,self.shiny); + obj = gluNewQuadric() + gluCylinder(obj,rad,rad,bond[10],self.nsides,self.nsides) + glPopMatrix() + + if self.tridraw: + fillflag = self.vizinfo.tfill[int(self.tridraw[0][1])] + + if fillflag != 1: + if fillflag: + glEnable(GL_POLYGON_OFFSET_FILL) + glPolygonOffset(1.0,1.0) + glBegin(GL_TRIANGLES) + ncolor = self.vizinfo.ntcolor + for tri in self.tridraw: + itype = int(tri[1]) + if itype > ncolor: raise StandardError,"tri type too big" + red,green,blue = self.vizinfo.tcolor[itype] + glMaterialfv(GL_FRONT_AND_BACK,GL_EMISSION,[red,green,blue,1.0]); + glMaterialf(GL_FRONT_AND_BACK,GL_SHININESS,self.shiny); + glNormal3f(tri[11],tri[12],tri[13]) + glVertex3f(tri[2],tri[3],tri[4]) + glVertex3f(tri[5],tri[6],tri[7]) + glVertex3f(tri[8],tri[9],tri[10]) + glEnd() + if fillflag: glDisable(GL_POLYGON_OFFSET_FILL) + + if fillflag: + glDisable(GL_LIGHTING) + glPolygonMode(GL_FRONT_AND_BACK,GL_LINE) + glLineWidth(self.bxthick) + glColor3f(self.bxcol[0],self.bxcol[1],self.bxcol[2]) + glBegin(GL_TRIANGLES) + for tri in self.tridraw: + glVertex3f(tri[2],tri[3],tri[4]) + glVertex3f(tri[5],tri[6],tri[7]) + glVertex3f(tri[8],tri[9],tri[10]) + glEnd() + glEnable(GL_LIGHTING) + glPolygonMode(GL_FRONT_AND_BACK,GL_FILL) + + # draw clipped scene = atoms, bonds, triangles + + else: + box = self.boxdraw + xlo = box[0] + self.clipxlo*(box[3] - box[0]) + xhi = box[0] + self.clipxhi*(box[3] - box[0]) + ylo = box[1] + self.clipylo*(box[4] - box[1]) + yhi = box[1] + self.clipyhi*(box[4] - box[1]) + zlo = box[2] + self.clipzlo*(box[5] - box[2]) + zhi = box[2] + self.clipzhi*(box[5] - box[2]) + + for atom in self.atomdraw: + x,y,z = atom[2],atom[3],atom[4] + if x >= xlo and x <= xhi and y >= ylo and y <= yhi and \ + z >= zlo and z <= zhi: + glTranslatef(x,y,z); + glCallList(self.calllist[int(atom[1])]); + glTranslatef(-x,-y,-z); + + if self.bonddraw: + bound = 0.25 * self.distance + ncolor = self.vizinfo.nbcolor + for bond in self.bonddraw: + xmin = min2(bond[2],bond[5]) + xmax = max2(bond[2],bond[5]) + ymin = min2(bond[3],bond[6]) + ymax = max2(bond[3],bond[6]) + zmin = min2(bond[4],bond[7]) + zmax = max2(bond[4],bond[7]) + if xmin >= xlo and xmax <= xhi and \ + ymin >= ylo and ymax <= yhi and zmin >= zlo and zmax <= zhi: + if bond[10] > bound: continue + itype = int(bond[1]) + if itype > ncolor: raise StandardError,"bond type too big" + red,green,blue = self.vizinfo.bcolor[itype] + rad = self.vizinfo.brad[itype] + glPushMatrix() + glTranslatef(bond[2],bond[3],bond[4]) + glRotatef(bond[11],bond[12],bond[13],0.0) + glMaterialfv(GL_FRONT_AND_BACK,GL_EMISSION,[red,green,blue,1.0]); + glMaterialf(GL_FRONT_AND_BACK,GL_SHININESS,self.shiny); + obj = gluNewQuadric() + gluCylinder(obj,rad,rad,bond[10],self.nsides,self.nsides) + glPopMatrix() + + if self.tridraw: + fillflag = self.vizinfo.tfill[int(self.tridraw[0][1])] + + if fillflag != 1: + if fillflag: + glEnable(GL_POLYGON_OFFSET_FILL) + glPolygonOffset(1.0,1.0) + glBegin(GL_TRIANGLES) + ncolor = self.vizinfo.ntcolor + for tri in self.tridraw: + xmin = min3(tri[2],tri[5],tri[8]) + xmax = max3(tri[2],tri[5],tri[8]) + ymin = min3(tri[3],tri[6],tri[9]) + ymax = max3(tri[3],tri[6],tri[9]) + zmin = min3(tri[4],tri[7],tri[10]) + zmax = max3(tri[4],tri[7],tri[10]) + if xmin >= xlo and xmax <= xhi and \ + ymin >= ylo and ymax <= yhi and \ + zmin >= zlo and zmax <= zhi: + itype = int(tri[1]) + if itype > ncolor: raise StandardError,"tri type too big" + red,green,blue = self.vizinfo.tcolor[itype] + glMaterialfv(GL_FRONT_AND_BACK,GL_EMISSION, + [red,green,blue,1.0]); + glMaterialf(GL_FRONT_AND_BACK,GL_SHININESS,self.shiny); + glNormal3f(tri[11],tri[12],tri[13]) + glVertex3f(tri[2],tri[3],tri[4]) + glVertex3f(tri[5],tri[6],tri[7]) + glVertex3f(tri[8],tri[9],tri[10]) + glEnd() + if fillflag: glDisable(GL_POLYGON_OFFSET_FILL) + + if fillflag: + glDisable(GL_LIGHTING) + glPolygonMode(GL_FRONT_AND_BACK,GL_LINE) + glLineWidth(self.bxthick) + glColor3f(self.bxcol[0],self.bxcol[1],self.bxcol[2]) + glBegin(GL_TRIANGLES) + for tri in self.tridraw: + xmin = min3(tri[2],tri[5],tri[8]) + xmax = max3(tri[2],tri[5],tri[8]) + ymin = min3(tri[3],tri[6],tri[9]) + ymax = max3(tri[3],tri[6],tri[9]) + zmin = min3(tri[4],tri[7],tri[10]) + zmax = max3(tri[4],tri[7],tri[10]) + if xmin >= xlo and xmax <= xhi and \ + ymin >= ylo and ymax <= yhi and \ + zmin >= zlo and zmax <= zhi: + glVertex3f(tri[2],tri[3],tri[4]) + glVertex3f(tri[5],tri[6],tri[7]) + glVertex3f(tri[8],tri[9],tri[10]) + glEnd() + glEnable(GL_LIGHTING) + glPolygonMode(GL_FRONT_AND_BACK,GL_FILL) + + if self.cache: glEndList() + + glFlush() + + # -------------------------------------------------------------------- + # make new call list for each atom type + # called when atom color/rad/quality is changed + + def make_atom_calllist(self): + # extend calllist array if necessary + + if self.vizinfo.nacolor > self.nclist: + for i in range(self.vizinfo.nacolor-self.nclist): self.calllist.append(0) + self.nclist = self.vizinfo.nacolor + + # create new calllist for each atom type + + for itype in xrange(1,self.vizinfo.nacolor+1): + if self.calllist[itype]: glDeleteLists(self.calllist[itype],1) + ilist = glGenLists(1) + self.calllist[itype] = ilist + glNewList(ilist,GL_COMPILE) + red,green,blue = self.vizinfo.acolor[itype] + rad = self.vizinfo.arad[itype] + glColor3f(red,green,blue); + +# glPointSize(10.0*rad) +# glBegin(GL_POINTS) +# glVertex3f(0.0,0.0,0.0) +# glEnd() + + glMaterialfv(GL_FRONT,GL_EMISSION,[red,green,blue,1.0]); + glMaterialf(GL_FRONT,GL_SHININESS,self.shiny); + glutSolidSphere(rad,self.nslices,self.nstacks) + glEndList() + + # -------------------------------------------------------------------- + # augment bond info returned by viz() with info needed for GL draw + # info = length, theta, -dy, dx for bond orientation + + def bonds_augment(self,bonds): + for bond in bonds: + dx = bond[5] - bond[2] + dy = bond[6] - bond[3] + dz = bond[7] - bond[4] + length = sqrt(dx*dx + dy*dy + dz*dz) + dx /= length + dy /= length + dz /= length + theta = acos(dz)*180.0/pi + bond += [length,theta,-dy,dx] + + # -------------------------------------------------------------------- + + def draw_box(self,flag): + xlo,ylo,zlo,xhi,yhi,zhi = self.boxdraw + + if flag: + tmp = xlo + self.clipxlo*(xhi - xlo) + xhi = xlo + self.clipxhi*(xhi - xlo) + xlo = tmp + tmp = ylo + self.clipylo*(yhi - ylo) + yhi = ylo + self.clipyhi*(yhi - ylo) + ylo = tmp + tmp = zlo + self.clipzlo*(zhi - zlo) + zhi = zlo + self.clipzhi*(zhi - zlo) + zlo = tmp + + glLineWidth(self.bxthick) + glColor3f(self.bxcol[0],self.bxcol[1],self.bxcol[2]) + + glBegin(GL_LINE_LOOP) + glVertex3f(xlo,ylo,zlo) + glVertex3f(xhi,ylo,zlo) + glVertex3f(xhi,yhi,zlo) + glVertex3f(xlo,yhi,zlo) + glEnd() + + glBegin(GL_LINE_LOOP) + glVertex3f(xlo,ylo,zhi) + glVertex3f(xhi,ylo,zhi) + glVertex3f(xhi,yhi,zhi) + glVertex3f(xlo,yhi,zhi) + glEnd() + + glBegin(GL_LINES) + glVertex3f(xlo,ylo,zlo) + glVertex3f(xlo,ylo,zhi) + glVertex3f(xhi,ylo,zlo) + glVertex3f(xhi,ylo,zhi) + glVertex3f(xhi,yhi,zlo) + glVertex3f(xhi,yhi,zhi) + glVertex3f(xlo,yhi,zlo) + glVertex3f(xlo,yhi,zhi) + glEnd() + + # -------------------------------------------------------------------- + + def draw_axes(self): + xlo,ylo,zlo,xhi,yhi,zhi = self.boxdraw + + delta = xhi-xlo + if yhi-ylo > delta: delta = yhi-ylo + if zhi-zlo > delta: delta = zhi-zlo + delta *= 0.1 + + glLineWidth(self.bxthick) + + glBegin(GL_LINES) + glColor3f(1,0,0) + glVertex3f(xlo-delta,ylo-delta,zlo-delta) + glVertex3f(xhi-delta,ylo-delta,zlo-delta) + glColor3f(0,1,0) + glVertex3f(xlo-delta,ylo-delta,zlo-delta) + glVertex3f(xlo-delta,yhi-delta,zlo-delta) + glColor3f(0,0,1) + glVertex3f(xlo-delta,ylo-delta,zlo-delta) + glVertex3f(xlo-delta,ylo-delta,zhi-delta) + glEnd() + + # -------------------------------------------------------------------- + + def save(self,file=None): + self.w.update() # force image on screen to be current before saving it + + pstring = glReadPixels(0,0,self.xpixels,self.ypixels, + GL_RGBA,GL_UNSIGNED_BYTE) + snapshot = Image.fromstring("RGBA",(self.xpixels,self.ypixels),pstring) + snapshot = snapshot.transpose(Image.FLIP_TOP_BOTTOM) + + if not file: file = self.file + snapshot.save(file + ".png") + + # -------------------------------------------------------------------- + + def adef(self): + self.vizinfo.setcolors("atom",range(100),"loop") + self.vizinfo.setradii("atom",range(100),0.45) + self.make_atom_calllist() + self.cachelist = -self.cachelist + self.w.tkRedraw() + + # -------------------------------------------------------------------- + + def bdef(self): + self.vizinfo.setcolors("bond",range(100),"loop") + self.vizinfo.setradii("bond",range(100),0.25) + self.cachelist = -self.cachelist + self.w.tkRedraw() + + # -------------------------------------------------------------------- + + def tdef(self): + self.vizinfo.setcolors("tri",range(100),"loop") + self.vizinfo.setfills("tri",range(100),0) + self.cachelist = -self.cachelist + self.w.tkRedraw() + + # -------------------------------------------------------------------- + + def ldef(self): + self.vizinfo.setcolors("line",range(100),"loop") + self.vizinfo.setradii("line",range(100),0.25) + self.cachelist = -self.cachelist + self.w.tkRedraw() + + # -------------------------------------------------------------------- + + def acol(self,atypes,colors): + self.vizinfo.setcolors("atom",atypes,colors) + self.make_atom_calllist() + self.cachelist = -self.cachelist + self.w.tkRedraw() + + # -------------------------------------------------------------------- + + def arad(self,atypes,radii): + self.vizinfo.setradii("atom",atypes,radii) + self.make_atom_calllist() + self.cachelist = -self.cachelist + self.w.tkRedraw() + + # -------------------------------------------------------------------- + + def bcol(self,btypes,colors): + self.vizinfo.setcolors("bond",btypes,colors) + self.cachelist = -self.cachelist + self.w.tkRedraw() + + # -------------------------------------------------------------------- + + def brad(self,btypes,radii): + self.vizinfo.setradii("bond",btypes,radii) + self.cachelist = -self.cachelist + self.w.tkRedraw() + + # -------------------------------------------------------------------- + + def tcol(self,ttypes,colors): + self.vizinfo.setcolors("tri",ttypes,colors) + self.cachelist = -self.cachelist + self.w.tkRedraw() + + # -------------------------------------------------------------------- + + def tfill(self,ttypes,flags): + self.vizinfo.setfills("tri",ttypes,flags) + self.cachelist = -self.cachelist + self.w.tkRedraw() + + # -------------------------------------------------------------------- + + def lcol(self,ltypes,colors): + self.vizinfo.setcolors("line",ltypes,colors) + self.cachelist = -self.cachelist + self.w.tkRedraw() + + # -------------------------------------------------------------------- + + def lrad(self,ltypes,radii): + self.vizinfo.setradii("line",ltypes,radii) + self.cachelist = -self.cachelist + self.w.tkRedraw() + +# -------------------------------------------------------------------- +# derived class from Togl's Opengl +# overwrite redraw, translate, rotate, scale methods +# latter 3 are mouse-motion methods + +class MyOpengl(Opengl): + def __init__(self, master, cnf={}, **kw): + args = (self,master,cnf) + Opengl.__init__(*args,**kw) + Opengl.autospin_allowed = 0 + + # redraw Opengl scene + # call parent redraw() method + + def tkRedraw(self,*dummy): + if not self.initialised: return + self.tk.call(self._w,'makecurrent') + self.redraw(self) + self.tk.call(self._w,'swapbuffers') + + # left button translate + # access parent xshift/yshift and call parent trans() method + + def tkTranslate(self,event): + dx = event.x - self.xmouse + dy = event.y - self.ymouse + x = self.parent.xshift + dx + y = self.parent.yshift - dy + self.parent.shift(x,y) + self.tkRedraw() + self.tkRecordMouse(event) + + # middle button trackball + # call parent mouse_rotate() method + + def tkRotate(self,event): + self.parent.mouse_rotate(event.x,event.y,self.xmouse,self.ymouse) + self.tkRedraw() + self.tkRecordMouse(event) + + # right button zoom + # access parent scale and call parent zoom() method + + def tkScale(self,event): + scale = 1 - 0.01 * (event.y - self.ymouse) + if scale < 0.001: scale = 0.001 + elif scale > 1000: scale = 1000 + scale *= self.parent.scale + self.parent.zoom(scale) + self.tkRedraw() + self.tkRecordMouse(event) + +# -------------------------------------------------------------------- +# draw a line segment + +def segment(p1,p2): + glVertex3f(p1[0],p1[1],p1[2]) + glVertex3f(p2[0],p2[1],p2[2]) + +# -------------------------------------------------------------------- +# normalize a 3-vector to unit length + +def vecnorm(v): + length = sqrt(v[0]*v[0] + v[1]*v[1] + v[2]*v[2]) + return [v[0]/length,v[1]/length,v[2]/length] + +# -------------------------------------------------------------------- +# dot product of two 3-vectors + +def vecdot(v1,v2): + return v1[0]*v2[0] + v1[1]*v2[1] + v1[2]*v2[2] + +# -------------------------------------------------------------------- +# cross product of two 3-vectors + +def veccross(v1,v2): + v = [0,0,0] + v[0] = v1[1]*v2[2] - v1[2]*v2[1] + v[1] = v1[2]*v2[0] - v1[0]*v2[2] + v[2] = v1[0]*v2[1] - v1[1]*v2[0] + return v + +# -------------------------------------------------------------------- +# return characteristic distance of simulation domain = max dimension + +def compute_distance(box): + distance = box[3]-box[0] + if box[4]-box[1] > distance: distance = box[4]-box[1] + if box[5]-box[2] > distance: distance = box[5]-box[2] + return distance + +# -------------------------------------------------------------------- +# return center of box as 3 vector + +def compute_center(box): + c = [0,0,0] + c[0] = 0.5 * (box[0] + box[3]) + c[1] = 0.5 * (box[1] + box[4]) + c[2] = 0.5 * (box[2] + box[5]) + return c + +# -------------------------------------------------------------------- +# return min of 2 values + +def min2(a,b): + if b < a: a = b + return a + +# -------------------------------------------------------------------- +# return max of 2 values + +def max2(a,b): + if b > a: a = b + return a + +# -------------------------------------------------------------------- +# return min of 3 values + +def min3(a,b,c): + if b < a: a = b + if c < a: a = c + return a + +# -------------------------------------------------------------------- +# return max of 3 values + +def max3(a,b,c): + if b > a: a = b + if c > a: a = c + return a diff --git a/python/examples/pizza/gnu.py b/python/examples/pizza/gnu.py new file mode 100644 index 0000000000..f6f0167330 --- /dev/null +++ b/python/examples/pizza/gnu.py @@ -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/NumPy 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 --git a/python/examples/pizza/vizinfo.py b/python/examples/pizza/vizinfo.py new file mode 100644 index 0000000000..acc421a1a6 --- /dev/null +++ b/python/examples/pizza/vizinfo.py @@ -0,0 +1,391 @@ +# 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. + +# vizinfo class, not a top-level Pizza.py tool + +# History +# 8/05, Matt Jones (BYU): original version +# 9/05, Steve Plimpton: added 140-color table + +# ToDo list + +# Variables + +# Imports and external programs + +import types + +# Class definition + +class vizinfo: + """ + Information holder for Pizza.py visualization tools + + acolor,bcolor,tcolor,lcolor = RGB values for each atom/bond/tri/line type + arad = radius of each atom type + brad,lrad = thickness of each bond/line type + tfill = fill flag for each triangle type + + all of these arrays are indexed by object type which runs 1-Ntype + nacolor,nbcolor,ntcolor,nlcolor,narad,nbrad,nlrad,ntfill + are # of types each array holds + actual length is nacolor+1 so that array can be indexed by 1-Ntype + + setcolors() = set atom/bond/tri/line colors + setradii() = set atom/bond/line radii/thickness + setfill() = set triangle fill factor + extend() = grow an array + """ + + # -------------------------------------------------------------------- + + def __init__(self): + self.acolor = [] + self.arad = [] + self.bcolor = [] + self.brad = [] + self.tcolor = [] + self.tfill = [] + self.lcolor = [] + self.lrad = [] + self.nacolor = self.narad = 0 + self.nbcolor = self.nbrad = 0 + self.ntcolor = self.ntfill = 0 + self.nlcolor = self.nlrad = 0 + + # -------------------------------------------------------------------- + # set color RGB for which = atoms, bonds, triangles + + def setcolors(self,which,ids,rgbs): + + # convert args into lists if single values + # if arg = 0, convert to full-range list + + if type(ids) is types.IntType and ids == 0: + if which == "atom": ids = range(self.nacolor) + if which == "bond": ids = range(self.nbcolor) + if which == "tri": ids = range(self.ntcolor) + if which == "line": ids = range(self.nlcolor) + if type(ids) is not types.ListType and type(ids) is not types.TupleType: + ids = [ids] + if type(rgbs) is not types.ListType and type(rgbs) is not types.TupleType: + rgbs = [rgbs] + + # if list of types has a 0, increment each type value + + if 0 in ids: + for i in xrange(len(ids)): ids[i] += 1 + + # extend storage list if necessary + # extend other arrays for same "which" so that gl::make_atom_calllist + # has valid arrays to work with + + if which == "atom": + if max(ids) > self.nacolor: + self.nacolor = self.extend(self.acolor,max(ids)) + self.nacolor = self.extend(self.arad,max(ids)) + if which == "bond": + if max(ids) > self.nbcolor: + self.nbcolor = self.extend(self.bcolor,max(ids)) + self.nbcolor = self.extend(self.brad,max(ids)) + if which == "tri": + if max(ids) > self.ntcolor: + self.ntcolor = self.extend(self.tcolor,max(ids)) + self.ntcolor = self.extend(self.tfill,max(ids)) + if which == "line": + if max(ids) > self.nlcolor: + self.nlcolor = self.extend(self.lcolor,max(ids)) + self.nlcolor = self.extend(self.lrad,max(ids)) + + # set color for each type + # if list lengths match, set directly, else interpolate + # convert final color from 0-255 to 0.0-1.0 + + ntypes = len(ids) + nrgbs = len(rgbs) + + for i in xrange(ntypes): + id = ids[i] + + if rgbs[0] == "loop": + list = colors.keys() + red,green,blue = colors[list[i % len(colors)]] + elif ntypes == nrgbs: + red,green,blue = colors[rgbs[i]] + else: + r = i/float(ntypes-1) * float(nrgbs-1) + jlo = int(r) + jhi = jlo + 1 + if jhi == nrgbs: jhi = nrgbs - 1 + clo = colors[rgbs[jlo]] + chi = colors[rgbs[jhi]] + delta = r - jlo + red = clo[0] + delta*(chi[0]-clo[0]) + green = clo[1] + delta*(chi[1]-clo[1]) + blue = clo[2] + delta*(chi[2]-clo[2]) + + color = [red/255.0,green/255.0,blue/255.0] + + if which == "atom": self.acolor[id] = color + if which == "bond": self.bcolor[id] = color + if which == "tri": self.tcolor[id] = color + if which == "line": self.lcolor[id] = color + + # -------------------------------------------------------------------- + # set radii for which = atoms, bonds, lines + + def setradii(self,which,ids,radii): + + # convert args into lists if single values + # if arg = 0, convert to full-range list + + if type(ids) is types.IntType and ids == 0: + if which == "atom": ids = range(self.narad) + if which == "bond": ids = range(self.nbrad) + if which == "line": ids = range(self.nlrad) + if type(ids) is not types.ListType and type(ids) is not types.TupleType: + ids = [ids] + if type(radii) is not types.ListType and \ + type(radii) is not types.TupleType: + radii = [radii] + + # if list of types has a 0, increment each type value + + if 0 in ids: + for i in xrange(len(ids)): ids[i] += 1 + + # extend storage list if necessary + # extend other arrays for same "which" so that gl::make_atom_calllist + # has valid arrays to work with + + if which == "atom": + if max(ids) > self.narad: + self.narad = self.extend(self.arad,max(ids)) + self.narad = self.extend(self.acolor,max(ids)) + if which == "bond": + if max(ids) > self.nbrad: + self.nbrad = self.extend(self.brad,max(ids)) + self.nbrad = self.extend(self.bcolor,max(ids)) + if which == "line": + if max(ids) > self.nlrad: + self.nlrad = self.extend(self.lrad,max(ids)) + self.nlrad = self.extend(self.lcolor,max(ids)) + + # set radius for each type + # if list lengths match, set directly, else interpolate + + ntypes = len(ids) + nradii = len(radii) + + for i in range(ntypes): + id = ids[i] + + if ntypes == nradii: rad = radii[i] + else: + r = i/float(ntypes-1) * float(nradii-1) + jlo = int(r) + jhi = jlo + 1 + if jhi == nradii: jhi = nradii - 1 + rlo = radii[jlo] + rhi = radii[jhi] + delta = r - jlo + rad = rlo + delta*(rhi-rlo) + + if which == "atom": self.arad[id] = rad + if which == "bond": self.brad[id] = rad + if which == "line": self.lrad[id] = rad + + # -------------------------------------------------------------------- + # set triangle fill style + # 0 = fill only, 1 = line only, 2 = fill and line + + def setfills(self,which,ids,fills): + + # convert args into lists if single values + # if arg = 0, convert to full-range list + + if type(ids) is types.IntType and ids == 0: + ids = range(self.ntfill) + if type(ids) is not types.ListType and type(ids) is not types.TupleType: + ids = [ids] + if type(fills) is not types.ListType and \ + type(fills) is not types.TupleType: + fills = [fills] + + # if list of types has a 0, increment each type value + + if 0 in ids: + for i in xrange(len(ids)): ids[i] += 1 + + # extend storage list if necessary + # extend other arrays for same "which" so that gl::make_atom_calllist + # has valid arrays to work with + + if max(ids) > self.ntfill: + self.ntfill = self.extend(self.tfill,max(ids)) + self.ntfill = self.extend(self.tcolor,max(ids)) + + # set fill flag for each type + # if list lengths match, set directly, else set types to 1st fill value + + if len(fills) == len(ids): + for i in xrange(len(ids)): self.tfill[ids[i]] = int(fills[i]) + else: + for id in ids: self.tfill[id] = int(fills[0]) + + # -------------------------------------------------------------------- + + def extend(self,array,n): + for i in range(n-len(array)+1): array.append(0) + return n + +# -------------------------------------------------------------------- +# dictionary of 140 color names and associated RGB values + +colors = {} + +colors["aliceblue"] = [240, 248, 255] +colors["antiquewhite"] = [250, 235, 215] +colors["aqua"] = [0, 255, 255] +colors["aquamarine"] = [127, 255, 212] +colors["azure"] = [240, 255, 255] +colors["beige"] = [245, 245, 220] +colors["bisque"] = [255, 228, 196] +colors["black"] = [0, 0, 0] +colors["blanchedalmond"] = [255, 255, 205] +colors["blue"] = [0, 0, 255] +colors["blueviolet"] = [138, 43, 226] +colors["brown"] = [165, 42, 42] +colors["burlywood"] = [222, 184, 135] +colors["cadetblue"] = [95, 158, 160] +colors["chartreuse"] = [127, 255, 0] +colors["chocolate"] = [210, 105, 30] +colors["coral"] = [255, 127, 80] +colors["cornflowerblue"] = [100, 149, 237] +colors["cornsilk"] = [255, 248, 220] +colors["crimson"] = [220, 20, 60] +colors["cyan"] = [0, 255, 255] +colors["darkblue"] = [0, 0, 139] +colors["darkcyan"] = [0, 139, 139] +colors["darkgoldenrod"] = [184, 134, 11] +colors["darkgray"] = [169, 169, 169] +colors["darkgreen"] = [0, 100, 0] +colors["darkkhaki"] = [189, 183, 107] +colors["darkmagenta"] = [139, 0, 139] +colors["darkolivegreen"] = [85, 107, 47] +colors["darkorange"] = [255, 140, 0] +colors["darkorchid"] = [153, 50, 204] +colors["darkred"] = [139, 0, 0] +colors["darksalmon"] = [233, 150, 122] +colors["darkseagreen"] = [143, 188, 143] +colors["darkslateblue"] = [72, 61, 139] +colors["darkslategray"] = [47, 79, 79] +colors["darkturquoise"] = [0, 206, 209] +colors["darkviolet"] = [148, 0, 211] +colors["deeppink"] = [255, 20, 147] +colors["deepskyblue"] = [0, 191, 255] +colors["dimgray"] = [105, 105, 105] +colors["dodgerblue"] = [30, 144, 255] +colors["firebrick"] = [178, 34, 34] +colors["floralwhite"] = [255, 250, 240] +colors["forestgreen"] = [34, 139, 34] +colors["fuchsia"] = [255, 0, 255] +colors["gainsboro"] = [220, 220, 220] +colors["ghostwhite"] = [248, 248, 255] +colors["gold"] = [255, 215, 0] +colors["goldenrod"] = [218, 165, 32] +colors["gray"] = [128, 128, 128] +colors["green"] = [0, 128, 0] +colors["greenyellow"] = [173, 255, 47] +colors["honeydew"] = [240, 255, 240] +colors["hotpink"] = [255, 105, 180] +colors["indianred"] = [205, 92, 92] +colors["indigo"] = [75, 0, 130] +colors["ivory"] = [255, 240, 240] +colors["khaki"] = [240, 230, 140] +colors["lavender"] = [230, 230, 250] +colors["lavenderblush"] = [255, 240, 245] +colors["lawngreen"] = [124, 252, 0] +colors["lemonchiffon"] = [255, 250, 205] +colors["lightblue"] = [173, 216, 230] +colors["lightcoral"] = [240, 128, 128] +colors["lightcyan"] = [224, 255, 255] +colors["lightgoldenrodyellow"] = [250, 250, 210] +colors["lightgreen"] = [144, 238, 144] +colors["lightgrey"] = [211, 211, 211] +colors["lightpink"] = [255, 182, 193] +colors["lightsalmon"] = [255, 160, 122] +colors["lightseagreen"] = [32, 178, 170] +colors["lightskyblue"] = [135, 206, 250] +colors["lightslategray"] = [119, 136, 153] +colors["lightsteelblue"] = [176, 196, 222] +colors["lightyellow"] = [255, 255, 224] +colors["lime"] = [0, 255, 0] +colors["limegreen"] = [50, 205, 50] +colors["linen"] = [250, 240, 230] +colors["magenta"] = [255, 0, 255] +colors["maroon"] = [128, 0, 0] +colors["mediumaquamarine"] = [102, 205, 170] +colors["mediumblue"] = [0, 0, 205] +colors["mediumorchid"] = [186, 85, 211] +colors["mediumpurple"] = [147, 112, 219] +colors["mediumseagreen"] = [60, 179, 113] +colors["mediumslateblue"] = [123, 104, 238] +colors["mediumspringgreen"] = [0, 250, 154] +colors["mediumturquoise"] = [72, 209, 204] +colors["mediumvioletred"] = [199, 21, 133] +colors["midnightblue"] = [25, 25, 112] +colors["mintcream"] = [245, 255, 250] +colors["mistyrose"] = [255, 228, 225] +colors["moccasin"] = [255, 228, 181] +colors["navajowhite"] = [255, 222, 173] +colors["navy"] = [0, 0, 128] +colors["oldlace"] = [253, 245, 230] +colors["olive"] = [128, 128, 0] +colors["olivedrab"] = [107, 142, 35] +colors["orange"] = [255, 165, 0] +colors["orangered"] = [255, 69, 0] +colors["orchid"] = [218, 112, 214] +colors["palegoldenrod"] = [238, 232, 170] +colors["palegreen"] = [152, 251, 152] +colors["paleturquoise"] = [175, 238, 238] +colors["palevioletred"] = [219, 112, 147] +colors["papayawhip"] = [255, 239, 213] +colors["peachpuff"] = [255, 239, 213] +colors["peru"] = [205, 133, 63] +colors["pink"] = [255, 192, 203] +colors["plum"] = [221, 160, 221] +colors["powderblue"] = [176, 224, 230] +colors["purple"] = [128, 0, 128] +colors["red"] = [255, 0, 0] +colors["rosybrown"] = [188, 143, 143] +colors["royalblue"] = [65, 105, 225] +colors["saddlebrown"] = [139, 69, 19] +colors["salmon"] = [250, 128, 114] +colors["sandybrown"] = [244, 164, 96] +colors["seagreen"] = [46, 139, 87] +colors["seashell"] = [255, 245, 238] +colors["sienna"] = [160, 82, 45] +colors["silver"] = [192, 192, 192] +colors["skyblue"] = [135, 206, 235] +colors["slateblue"] = [106, 90, 205] +colors["slategray"] = [112, 128, 144] +colors["snow"] = [255, 250, 250] +colors["springgreen"] = [0, 255, 127] +colors["steelblue"] = [70, 130, 180] +colors["tan"] = [210, 180, 140] +colors["teal"] = [0, 128, 128] +colors["thistle"] = [216, 191, 216] +colors["tomato"] = [253, 99, 71] +colors["turquoise"] = [64, 224, 208] +colors["violet"] = [238, 130, 238] +colors["wheat"] = [245, 222, 179] +colors["white"] = [255, 255, 255] +colors["whitesmoke"] = [245, 245, 245] +colors["yellow"] = [255, 255, 0] +colors["yellowgreen"] = [154, 205, 50] diff --git a/python/examples/plot.py b/python/examples/plot.py new file mode 100755 index 0000000000..6fd0636443 --- /dev/null +++ b/python/examples/plot.py @@ -0,0 +1,75 @@ +#!/usr/local/bin/python -i +# preceeding line should have path for Python on your machine + +# plot.py +# Purpose: plot Temp of running LAMMPS simulation via GnuPlot in Pizza.py +# Syntax: plot.py in.lammps Nfreq Nsteps compute-ID +# in.lammps = LAMMPS input script +# Nfreq = plot data point every this many steps +# Nsteps = run for this many steps +# compute-ID = ID of compute that calculates temperature +# (or any other scalar quantity) + +import sys +sys.path.append("./pizza") +from gnu import gnu + +# parse command line + +argv = sys.argv +if len(argv) != 5: + print "Syntax: plot.py in.lammps Nfreq Nsteps compute-ID" + sys.exit() + +infile = sys.argv[1] +nfreq = int(sys.argv[2]) +nsteps = int(sys.argv[3]) +compute = sys.argv[4] + +me = 0 +# uncomment if running in parallel via Pypar +#import pypar +#me = pypar.rank() +#nprocs = pypar.size() + +from lammps import lammps +lmp = lammps() + +# run infile all at once +# assumed to have no run command in it + +lmp.file(infile) +lmp.command("thermo %d" % nfreq) + +# initial 0-step run to generate initial 1-point plot + +lmp.command("run 0 pre yes post no") +value = lmp.extract_compute(compute,0,0) +ntimestep = 0 +xaxis = [ntimestep] +yaxis = [value] + +# wrapper on GnuPlot via Pizza.py gnu tool +# just proc 0 handles plotting + +if me == 0: + gn = gnu() + gn.plot(xaxis,yaxis) + gn.xrange(0,nsteps) + gn.title(compute,"Timestep","Temperature") + +# run nfreq steps at a time w/out pre/post, query compute, refresh plot + +while ntimestep < nsteps: + lmp.command("run %d pre no post no" % nfreq) + ntimestep += nfreq + value = lmp.extract_compute(compute,0,0) + xaxis.append(ntimestep) + yaxis.append(value) + if me == 0: gn.plot(xaxis,yaxis) + +lmp.command("run 0 pre no post yes") + +# uncomment if running in parallel via Pypar +#print "Proc %d out of %d procs has" % (me,nprocs), lmp +#pypar.finalize() diff --git a/python/examples/simple.py b/python/examples/simple.py new file mode 100755 index 0000000000..3505f9c707 --- /dev/null +++ b/python/examples/simple.py @@ -0,0 +1,50 @@ +#!/usr/local/bin/python -i +# preceeding line should have path for Python on your machine + +# simple.py +# Purpose: mimic operation of couple/simple/simple.cpp via Python +# Syntax: simple.py in.lammps +# in.lammps = LAMMPS input script + +import sys + +# parse command line + +argv = sys.argv +if len(argv) != 2: + print "Syntax: simple.py in.lammps" + sys.exit() + +infile = sys.argv[1] + +me = 0 +# uncomment if running in parallel via Pypar +#import pypar +#me = pypar.rank() +#nprocs = pypar.size() + +from lammps import lammps +lmp = lammps() + +# run infile one line at a time + +lines = open(infile,'r').readlines() +for line in lines: lmp.command(line) + +# run 10 more steps +# get coords from LAMMPS +# change coords of 1st atom +# put coords back into LAMMPS +# run a single step with changed coords + +lmp.command("run 10") +x = lmp.get_coords() +epsilon = 0.1 +x[0] += epsilon +lmp.put_coords(x) +lmp.command("run 1"); +lmp.command("run 1") + +# uncomment if running in parallel via Pypar +#print "Proc %d out of %d procs has" % (me,nprocs), lmp +#pypar.finalize() diff --git a/python/examples/trivial.py b/python/examples/trivial.py new file mode 100755 index 0000000000..644e2cbf43 --- /dev/null +++ b/python/examples/trivial.py @@ -0,0 +1,40 @@ +#!/usr/local/bin/python -i +# preceeding line should have path for Python on your machine + +# trivial.py +# Purpose: run a LAMMPS input script via Python +# Syntax: trivial.py in.lammps +# in.lammps = LAMMPS input script + +import sys + +# parse command line + +argv = sys.argv +if len(argv) != 2: + print "Syntax: trivial.py in.lammps" + sys.exit() + +infile = sys.argv[1] + +me = 0 +# uncomment if running in parallel via Pypar +#import pypar +#me = pypar.rank() +#nprocs = pypar.size() + +from lammps import lammps +lmp = lammps() + +# run infile all at once + +lmp.file(infile) + +# run infile one line at a time + +#lines = open(infile,'r').readlines() +#for line in lines: lmp.command(line) + +# uncomment if running in parallel via Pypar +#print "Proc %d out of %d procs has" % (me,nprocs), lmp +#pypar.finalize() diff --git a/python/examples/viz.py b/python/examples/viz.py new file mode 100755 index 0000000000..16b056123e --- /dev/null +++ b/python/examples/viz.py @@ -0,0 +1,82 @@ +#!/usr/local/bin/python -i +# preceeding line should have path for Python on your machine + +# viz.py +# Purpose: viz running LAMMPS simulation via GL tool in Pizza.py +# Syntax: viz.py in.lammps Nfreq Nsteps +# in.lammps = LAMMPS input script +# Nfreq = dump and viz shapshot every this many steps +# Nsteps = run for this many steps + +import sys +sys.path.append("./pizza") + +# parse command line + +argv = sys.argv +if len(argv) != 4: + print "Syntax: viz.py in.lammps Nfreq Nsteps" + sys.exit() + +infile = sys.argv[1] +nfreq = int(sys.argv[2]) +nsteps = int(sys.argv[3]) + +me = 0 +# uncomment if running in parallel via Pypar +#import pypar +#me = pypar.rank() +#nprocs = pypar.size() + +from lammps import lammps +lmp = lammps() + +# run infile all at once +# assumed to have no run command in it + +lmp.file(infile) +lmp.command("thermo %d" % nfreq) +lmp.command("dump python all atom %d tmp.dump" % nfreq) + +# initial 0-step run to generate dump file and image + +lmp.command("run 0 pre yes post no") +ntimestep = 0 + +# wrapper on GL window via Pizza.py gl tool +# just proc 0 handles reading of dump file and viz + +if me == 0: + import Tkinter + tkroot = Tkinter.Tk() + tkroot.withdraw() + + from dump import dump + from gl import gl + + d = dump("tmp.dump",0) + g = gl(d) + d.next() + d.unscale() + g.zoom(1) + g.shift(0,0) + g.rotate(0,270) + g.q(10) + g.box(1) + g.show(ntimestep) + +# run nfreq steps at a time w/out pre/post, read dump snapshot, display it + +while ntimestep < nsteps: + lmp.command("run %d pre no post no" % nfreq) + ntimestep += nfreq + if me == 0: + d.next() + d.unscale() + g.show(ntimestep) + +lmp.command("run 0 pre no post yes") + +# uncomment if running in parallel via Pypar +#print "Proc %d out of %d procs has" % (me,nprocs), lmp +#pypar.finalize() diff --git a/python/examples/vizplotgui.py b/python/examples/vizplotgui.py new file mode 100755 index 0000000000..eb436e3cb7 --- /dev/null +++ b/python/examples/vizplotgui.py @@ -0,0 +1,171 @@ +#!/usr/local/bin/python -i +# preceeding line should have path for Python on your machine + +# vizplotgui.py +# Purpose: viz running LAMMPS simulation with plot and GUI +# Syntax: vizplotgui.py in.lammps Nfreq compute-ID +# in.lammps = LAMMPS input script +# Nfreq = plot data point and viz shapshot every this many steps +# compute-ID = ID of compute that calculates temperature +# (or any other scalar quantity) + +import sys,time +sys.path.append("./pizza") + +# methods called by GUI + +def run(): + global runflag + runflag = 1 +def stop(): + global runflag + runflag = 0 +def settemp(value): + global temptarget + temptarget = slider.get() +def quit(): + global breakflag + breakflag = 1 + +# method called by timestep loop every Nfreq steps +# read dump snapshot and viz it, update plot with compute value + +def update(ntimestep): + d.next() + d.unscale() + g.show(ntimestep) + value = lmp.extract_compute(compute,0,0) + xaxis.append(ntimestep) + yaxis.append(value) + gn.plot(xaxis,yaxis) + +# parse command line + +argv = sys.argv +if len(argv) != 4: + print "Syntax: vizplotgui.py in.lammps Nfreq compute-ID" + sys.exit() + +infile = sys.argv[1] +nfreq = int(sys.argv[2]) +compute = sys.argv[3] + +me = 0 +# uncomment if running in parallel via Pypar +#import pypar +#me = pypar.rank() +#nprocs = pypar.size() + +from lammps import lammps +lmp = lammps() + +# run infile all at once +# assumed to have no run command in it + +lmp.file(infile) +lmp.command("thermo %d" % nfreq) +lmp.command("dump python all atom %d tmp.dump" % nfreq) + +# initial 0-step run to generate initial 1-point plot, dump file, and image + +lmp.command("run 0 pre yes post no") +value = lmp.extract_compute(compute,0,0) +ntimestep = 0 +xaxis = [ntimestep] +yaxis = [value] + +# wrapper on GL window via Pizza.py gl tool +# just proc 0 handles reading of dump file and viz + +breakflag = 0 +runflag = 0 +temptarget = 1.0 + +if me == 0: + from Tkinter import * + tkroot = Tk() + tkroot.withdraw() + + from dump import dump + from gl import gl + + d = dump("tmp.dump",0) + g = gl(d) + d.next() + d.unscale() + g.zoom(1) + g.shift(0,0) + g.rotate(0,270) + g.q(10) + g.box(1) + g.show(ntimestep) + +# display GUI with run/stop buttons and slider for temperature + +if me == 0: + from Tkinter import * + tkroot = Tk() + tkroot.withdraw() + root = Toplevel(tkroot) + root.title("LAMMPS GUI") + + frame = Frame(root) + Button(frame,text="Run",command=run).pack(side=LEFT) + Button(frame,text="Stop",command=stop).pack(side=LEFT) + slider = Scale(frame,from_=0.0,to=5.0,resolution=0.1, + orient=HORIZONTAL,label="Temperature") + slider.bind('',settemp) + slider.set(temptarget) + slider.pack(side=LEFT) + Button(frame,text="Quit",command=quit).pack(side=RIGHT) + frame.pack() + tkroot.update() + +# wrapper on GnuPlot via Pizza.py gnu tool + +if me == 0: + from gnu import gnu + gn = gnu() + gn.plot(xaxis,yaxis) + gn.title(compute,"Timestep","Temperature") + +# endless loop, checking status of GUI settings every Nfreq steps +# run with pre yes/no and post yes/no depending on go/stop status +# re-invoke fix langevin with new seed when temperature slider changes +# after re-invoke of fix langevin, run with pre yes + +running = 0 +temp = temptarget +seed = 12345 + +lmp.command("fix 2 all langevin %g %g 0.1 %d" % (temp,temp,seed)) + +while 1: + if me == 0: tkroot.update() + if temp != temptarget: + temp = temptarget + seed += me+1 + lmp.command("fix 2 all langevin %g %g 0.1 12345" % (temp,temp)) + running = 0 + if runflag and running: + lmp.command("run %d pre no post no" % nfreq) + ntimestep += nfreq + if me == 0: update(ntimestep) + elif runflag and not running: + lmp.command("run %d pre yes post no" % nfreq) + ntimestep += nfreq + if me == 0: update(ntimestep) + elif not runflag and running: + lmp.command("run %d pre no post yes" % nfreq) + ntimestep += nfreq + if me == 0: update(ntimestep) + if breakflag: break + if runflag: running = 1 + else: running = 0 + time.sleep(0.01) + +lmp.command("run 0 pre no post yes") + +# uncomment if running in parallel via Pypar +#print "Proc %d out of %d procs has" % (me,nprocs), lmp +#pypar.finalize() diff --git a/python/lammps.py b/python/lammps.py new file mode 100644 index 0000000000..5f5ee52636 --- /dev/null +++ b/python/lammps.py @@ -0,0 +1,169 @@ +# ---------------------------------------------------------------------- +# LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator +# http://lammps.sandia.gov, Sandia National Laboratories +# Steve Plimpton, sjplimp@sandia.gov +# +# Copyright (2003) 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. +# +# See the README file in the top-level LAMMPS directory. +# ------------------------------------------------------------------------- + +# Python wrapper on LAMMPS library via ctypes + +import types +from ctypes import * + +LMPINT = 0 +LMPDOUBLE = 1 +LMPIPTR = 2 +LMPDPTR = 3 +LMPDPTRPTR = 4 + +class lammps: + def __init__(self,args=None): + + # attempt to load parallel library first, serial library next + # could provide caller a flag to choose which library to load + + try: + self.lib = CDLL("_lammps.so") + except: + try: + self.lib = CDLL("_lammps_serial.so") + except: + raise StandardError,"Could not load LAMMPS dynamic library" + + # create an instance of LAMMPS + # don't know how to pass an MPI communicator from PyPar + # no_mpi call lets LAMMPS use MPI_COMM_WORLD + # cargs = array of C strings from args + + if args: + args.insert(0,"lammps.py") + narg = len(args) + cargs = (c_char_p*narg)(*args) + self.lmp = c_void_p() + self.lib.lammps_open_no_mpi(narg,cargs,byref(self.lmp)) + else: + self.lmp = c_void_p() + self.lib.lammps_open_no_mpi(0,None,byref(self.lmp)) + # could use just this if LAMMPS lib interface supported it + # self.lmp = self.lib.lammps_open_no_mpi(0,None) + + def __del__(self): + if self.lmp: self.lib.lammps_close(self.lmp) + + def close(self): + self.lib.lammps_close(self.lmp) + self.lmp = None + + def file(self,file): + self.lib.lammps_file(self.lmp,file) + + def command(self,cmd): + self.lib.lammps_command(self.lmp,cmd) + + def extract_global(self,name,type): + if type == LMPDOUBLE: + self.lib.lammps_extract_global.restype = POINTER(c_double) + ptr = self.lib.lammps_extract_global(self.lmp,name) + return ptr[0] + if type == LMPINT: + self.lib.lammps_extract_global.restype = POINTER(c_int) + ptr = self.lib.lammps_extract_global(self.lmp,name) + return ptr[0] + return None + + def extract_atom(self,name,type): + if type == LMPDPTRPTR: + self.lib.lammps_extract_atom.restype = POINTER(POINTER(c_double)) + ptr = self.lib.lammps_extract_atom(self.lmp,name) + return ptr + if type == LMPDPTR: + self.lib.lammps_extract_atom.restype = POINTER(c_double) + ptr = self.lib.lammps_extract_atom(self.lmp,name) + return ptr + if type == LMPIPTR: + self.lib.lammps_extract_atom.restype = POINTER(c_int) + ptr = self.lib.lammps_extract_atom(self.lmp,name) + return ptr + return None + + def extract_compute(self,id,style,type): + if type == 0: + if style > 0: return None + self.lib.lammps_extract_compute.restype = POINTER(c_double) + ptr = self.lib.lammps_extract_compute(self.lmp,id,style,type) + return ptr[0] + elif type == 1: + self.lib.lammps_extract_compute.restype = POINTER(c_double) + ptr = self.lib.lammps_extract_compute(self.lmp,id,style,type) + return ptr + elif type == 2: + self.lib.lammps_extract_compute.restype = POINTER(POINTER(c_double)) + ptr = self.lib.lammps_extract_compute(self.lmp,id,style,type) + return ptr + return None + + # in case of global datum, free memory for 1 double via lammps_free() + # double was allocated by library interface function + + def extract_fix(self,id,style,type,i=0,j=0): + if type == 0: + if style > 0: return None + self.lib.lammps_extract_fix.restype = POINTER(c_double) + ptr = self.lib.lammps_extract_bix(self.lmp,id,style,type,i,j) + result = ptr[0] + self.lib.lammps_free(ptr) + return result + elif type == 1: + self.lib.lammps_extract_fix.restype = POINTER(c_double) + ptr = self.lib.lammps_extract_fix(self.lmp,id,style,type,i,j) + return ptr + elif type == 2: + self.lib.lammps_extract_fix.restype = POINTER(POINTER(c_double)) + ptr = self.lib.lammps_extract_fix(self.lmp,id,style,type,i,j) + return ptr + return None + + # free memory for 1 double or 1 vector of doubles via lammps_free() + # for vector, must copy nlocal returned values to local c_double vector + # memory was allocated by library interface function + + def extract_variable(self,name,group,type): + if type == 0: + self.lib.lammps_extract_variable.restype = POINTER(c_double) + ptr = self.lib.lammps_extract_variable(self.lmp,name,group) + result = ptr[0] + self.lib.lammps_free(ptr) + return result + if type == 1: + self.lib.lammps_extract_global.restype = POINTER(c_int) + nlocalptr = self.lib.lammps_extract_global(self.lmp,"nlocal") + nlocal = nlocalptr[0] + result = (c_double*nlocal)() + self.lib.lammps_extract_variable.restype = POINTER(c_double) + ptr = self.lib.lammps_extract_variable(self.lmp,name,group) + for i in xrange(nlocal): result[i] = ptr[i] + self.lib.lammps_free(ptr) + return result + return None + + def get_natoms(self): + return self.lib.lammps_get_natoms(self.lmp) + + def get_coords(self): + nlen = 3 * self.lib.lammps_get_natoms(self.lmp) + coords = (c_double*nlen)() + self.lib.lammps_get_coords(self.lmp,coords) + return coords + + # assume coords is an array of c_double, as created by get_coords() + # could check if it is some other Python object and create c_double array? + # constructor for c_double array can take an arg to use to fill it? + + def put_coords(self,coords): + self.lib.lammps_put_coords(self.lmp,coords) diff --git a/python/setup.py b/python/setup.py new file mode 100755 index 0000000000..80fcfcdb20 --- /dev/null +++ b/python/setup.py @@ -0,0 +1,39 @@ +#!/usr/local/bin/python + +""" +setup.py file for LAMMPS with system MPI library +""" + +from distutils.core import setup, Extension + +import os, glob +path = os.path.dirname(os.getcwd()) + +# list of src files for LAMMPS + +libfiles = glob.glob("%s/src/*.cpp" % path) + +lammps_library = Extension("_lammps", + sources = libfiles, + define_macros = [("MPICH_IGNORE_CXX_SEEK",1), + ("LAMMPS_GZIP",1), + ("FFT_NONE",1),], + # src files for LAMMPS + include_dirs = ["../src"], + # additional libs for MPICH on Linux + libraries = ["mpich","rt"], + # where to find the MPICH lib on Linux + library_dirs = ["/usr/local/lib"], + # additional libs for MPI on Mac + # libraries = ["mpi"], + ) + +setup(name = "lammps", + version = "26Oct10", + author = "Steve Plimpton", + author_email = "sjplimp@sandia.gov", + url = "http://lammps.sandia.gov", + description = """LAMMPS molecular dynamics library - parallel""", + py_modules = ["lammps"], + ext_modules = [lammps_library] + ) diff --git a/python/setup_serial.py b/python/setup_serial.py new file mode 100755 index 0000000000..425f47062b --- /dev/null +++ b/python/setup_serial.py @@ -0,0 +1,34 @@ +#!/usr/local/bin/python + +""" +setup_serial.py file for LAMMPS with dummy serial MPI library +""" + +from distutils.core import setup, Extension + +import os, glob +path = os.path.dirname(os.getcwd()) + +# list of src files for LAMMPS and MPI STUBS + +libfiles = glob.glob("%s/src/*.cpp" % path) + \ + glob.glob("%s/src/STUBS/*.cpp" % path) + +lammps_library = Extension("_lammps_serial", + sources = libfiles, + define_macros = [("MPICH_IGNORE_CXX_SEEK",1), + ("LAMMPS_GZIP",1), + ("FFT_NONE",1),], + # src files for LAMMPS and MPI STUBS + include_dirs = ["../src", "../src/STUBS"] + ) + +setup(name = "lammps_serial", + version = "26Oct10", + author = "Steve Plimpton", + author_email = "sjplimp@sandia.gov", + url = "http://lammps.sandia.gov", + description = """LAMMPS molecular dynamics library - serial""", + py_modules = ["lammps"], + ext_modules = [lammps_library] + )